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PREFACE 


Bounded  nominal  paths  can  be  constructed  in  the  vicinity  of  the 
interior  equilibrium  point  (sometimes  called  a  libration  or  Lagrange 
point)  for  the  Sun-Ear th+Moon  Elliptic  Restricted  Three-Body  Problem. 
Numerical  integration  is  used  to  generate  the  periodic  or  quasi- 
per iodic  reference  trajectories  in  this  effort,  and  the  numerical  data 
is  then  curve  fit  using  a  cubic  spline  routine.  The  force  model  used 
in  this  effort  includes  the  gravitational  attractions  of  the  Sun  plus 
the  Ear th+Moon  barycenter  (treated  as  one  body),  along  with  the 
centrifugal  force  associated  with  rotation  of  the  system.  Solar 
radiation  pressure  force  is  also  included.  A  spacecraft  in  a  libration 
point  orbit  between  the  Earth  and  the  Sun  can  study  the  interaction  of 
the  Sun’s  corona  with  the  terrestrial  environment  and  will  be  of  great 
scientific  value. 

However,  the  spacecraft  in  such  an  orbit  will  drift  from  the 
nominal  (unstable)  path,  and  the  forces  individually  have  some  level  of 
uncertainty.  Both  range  and  range-rate  tracking  also  include  some 
inaccuracy  in  the  data  obtained.  The  accumulated  error  in  the 
spacecraft’s  position  and  velocity  relative  to  the  nominal  path  after  a 
predetermined  period  of  tracking  can  be  measured.  This  error,  or 
uncertainty,  in  the  spacecraft  state  is  measured  through  simulations 
commonly  referred  to  as  orbit  determination  error  analysis.  The 
specific  error  analysis  method  used  here  is  covariance  analysis;  the 
covariance  matrix,  after  a  predetermined  period  of  tracking  along  the 
nominal  path,  contains  the  state  uncertainties  (state  vector  element 
variances)  along  its  diagonal. 

This  effort  is  supported  by  the  Frank  J.  Seiler  Research 
Laboratory  and  has  been  conducted  as  doctoral  research  under  the 
direction  of  Professor  K.C.  Howell,  School  of  Aeronautics  and 
Astronautics,  Purdue  University,  West  Lafayette,  Indiana. 
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INTRODUCTION 


With  the  expansion  of  space  exploration  programs  worldwide, 
interest  has  incre;  jed  in  the  design  of  innovative,  complex,  and  yet 
low-cost  spacecraft  trajectories  that  meet  demanding  mission 
requirements.  In  most  of  the  missions  flown  in  the  last  few  decades, 
the  spacecraft  spent  the  majority  of  the  flight  time  in  a  force 
environment  dominated  by  a  single  gravitational  field.  For  the 
preliminary  mission  analysis  in  these  cases,  additional  attracting 
bodies  and  other  forces  could  be  modeled,  when  required,  as  perturbing 
influences.  An?'’ysis  of  some  recently  proposed  and  more  adventurous 
missions,  such  as  those  involving  libratlon  point  orbits,  will  require 
dynamic  r<ddels  of  higher  complexity,  since  at  least  two  gravitational 
fields  are'  of  nearly  equal  influence  on  the  spacecraft  throughout  the 
majority  of  the  mission.  Thus,  trajectories  determined  for  a  system 
consisting  of  numerous  gravitational  forces  have  been  of  particular 
theoretical  and  practical  interest  in  recent  years. 

One  type  of  many-body  problem,  motion  within  a  three-body  system 
of  forces,  has  a  wide  range  of  applications.  The  general  problem  of 
three  bodies  assumes  that  each  body  has  finite  mass  and  that  the  motion 
is  a  result  of  mutual  gravitational  attraction.  When  the  mass  of  one 
of  the  three  bodies  is  assumed  to  be  sufficiently  small  (infinitesimal) 
so  that  it  does  not  affect  the  motion  of  the  other  two  bodies 
(primaries)  in  the  system,  the  "restricted  three-body  problem"  results. 
The  primaries  can  be  further  assumed  to  be  moving  in  known  elliptic  or 
circular  orbits  about  their  common  center  of  mass.  Therefore,  the 
elliptic  restricted  three-body  problem,  where  the  primaries  are  assumed 
to  be  in  known  elliptic  orbits,  may  be  considered  a  reasonably 
approximate  model  for  a  spacecraft  movi  ^  within  the  gravitational 
fields  of  the  Sun  and  the  Earth,  for  instance. 
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In  the  formulation  of  the  restricted  three-body  problem,  one  mass 
is  defined  as  infinitesimal  relative  to  the  remaining  two  masses 
(primaries).  The  primaries,  unaffected  by  the  infinitesimal  mass,  move 
under  their  mutual  gravitational  attractions.  In  the  elliptic 
restricted  three-body  problem  (ER3BP),  the  primaries  are  assumed  to 
move  on  elliptic  paths.  If  the  eccentricity  of  the  primaries'  orbit  is 
assumed  to  be  zero,  the  circular  restricted  three-body  problem  (CR3BP) 
results.  Even  for  known  primary  motion,  however,  a  general, 
closed-form  solution  for  motion  of  the  third  body  of  infinitesimal  mass 
does  not  exist.  In  the  restricted  three-body  problem  (ro3BP  or  CE3BP), 
five  equilibrium  (libration)  solutions  can  be  found.  These  equilibrium 
points,  sometimes  also  referred  to  as  Lagrange  points,  are  particular 
solutions  of  the  equations  of  motion  governing  the  path  of  the 
infinitesimal  mass  moving  within  the  gravitational  fields  of  the 
primaries. 

The  equilibrium  points  are  defined  relative  to  a  coordinate  system 
rotating  with  the  primaries.  At  these  locations,  the  forces  on  the 
spacecraft  are  in  equilibrium.  These  forces  include  the  gravitational 
forces  from  the  massive  bodies  and  the  centrifugal  force  associated 
with  the  rotation  of  the,  system.  (The  addition  of  solar  radiation 
pressure  to  the  force  model  changes  the  locations  of  the  five  Lagrange 
points,  although  they  can  still  be  defined,  and  these  solar  radiation 
effects  are  discussed  in  Gordon^*^.)  The  libration  points  are  located 
in  the  plane  of  primary  rotation.  Three  of  the  libration  points  are  on 
the  line  between  the  two  massive  bodies,  and  one  of  these  col  linear 
points  is  interior  to  the  primaries.  The  last  two  points  are  at  the 
vertices  of  two  equilateral  triangles  in  the  plane  of  primary  rotation. 
The  triangles  have  a  common  base  that  is  the  line  between  the  primary 
masses. 

For  the  CR3BP,  the  five  libration  points  are  stationary  relative 
to  the  rotating  reference  frame.  If  the  problem  is  generalized  to  the 
ER3BP,  the  libration  points  pulsate  as  the  distance  between  the 
primaries  varies  with  time.  In  both  the  circular  and  elliptic 
restricted  problems,  two-dimensional  and  three-dimensional 
trajectories,  both  periodic  and  quasi-periodic  paths,  can  be  compiuted 
in  the  vicinity  of  these  libration  points. 
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Three-dimensional,  periodic  “halo"  orbits  in  the  vicinity  of  the 

collinear  libration  points  have  been  studied  since  the  late  1960s. 

Early  work  concerning  these  orbits  was  motivated  by  studies  related  to 

exploring  the  far  side  of  the  Moon.  These  studies  were  completed  in 

support  of  the  planned  Apollo  18  lunar  exploration  mission  that  was 

later  canceled.  Robert  Farquhar  coined  the  term  "halo"  to  describe  a 

three-dimensional,  periodic  orbit  near  a  libration  point  on  the  far 

121 

side  of  "the  Moon  in  the  Earth-Moon  system.  These  orbits  were 

designed- to  be  large  enough  so  that  the  spacecraft  would  be  cons’^antly 
in  view  of  the  Earth  and  would  thus  appear  as  a  halo  around  the  Moon. 
A  communications  station  in  this  type  of  orbit  could  maintain  constant 
contact  between  the  Earth  and  a  lunar  experimentation  station  on  the 
far  side  of  the  Moon. 

Quasi-periodic  orbits  near  libration  points  are  also  currently  of 

great  research  interest.  The  variations  in  size  and  shape  that  a 

quasi-periodic  orbit  can  exhibit  may  add  valuable  flexibility  for 

mission  planning.  This  type  of  bounded,  three-dimensional  libration 

point  trajectory  is  called  a  Lissajous  orbit  since  specific  planar 

projections  of  these  quasi-periodic  trajectories  may  look  like  a 

special  type  of  "Lissajous"  curve.  Physicist  Jules  Antoine  Lissajous 

(1822-1880)  investigated  curves  that  were  generated  by  compounding 

simple  harmonic  motions  at  right  angles,  and  he  delivered  a  paper  on 

this  subject  to  the  Paris  Academy  of  Sciences  in  1857.  Nathaniel 

Bowditch  of  Salem,  Massachusetts,  had  conducted  some  similar  work  in 

1815.  Lissajous  curves  have  a  wide  variety  of  shapes  that  depend  on 

the  frequency,  phase,  and  amplitude  of  the  orthogonal  components  of  the 
(4  5] 

motion.  ’  When  the  in-plane  and  the  (orthogonal)  out-of-plane 
frequencies  of  the  linearized  motion  are  nearly  (but  not)  equal,  the 
resulting  path  is  typically  called  a  Lissajous  trajectory. 

A  method  to  generate  approximations  for  this  type  of 
quasi-periodic  orbital  path  was  developed  analytically  by  Farquhar  and 
Kamel  in  1973.*^^  They  derived  a  third-order  approximate  analytic 
solution  for  a  translunar  libration  point  orbit  in  the  Earth-Moon  ER3BP 
that  also  included  solar  gravity  perturbations.  In  1975,  Richardson 
and  Cary  then  developed  a  fourth-order  analytic  Lissajous  approximation 
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[71 

in  the  Sun-Earth+Moon  barycenter  system.  The  notation  "Earth+Moon" 

indicates  that  the  Earth  and  the  Moon  are  treated  as  one  body  with  mass 

center  at  the  Earth-Moon  barycenter.  In  consideration  of  the  lunar 

perturbation,  Farquhar  has  shown  that  the  accuracy  of  solutions  in  the 

Sun-Earth  CR3BP  can  be  enhanced  if  the  collinear  libration  points  are 

defined  along  the  line  between  the  Sun  and  the  center  of  mass  of  the 

Earth  and  the  Moon.  Since  1975,  a  few  researchers  have  considered 

methods  to  numerically  generate  Lissajous  trajectories,  but  the  lack  of 

periodicity  of  a  Lissajous  path  complicates  numerical  construction  of 

bounded  trajectories.  Howell  and  Pernicka  have  developed  a  numerical 

technique  for  determination  of  three-dimensional,  bounded  Lissajous 

[9-t4] 

trajectories  of  arbitrary  size  and  duration.  Orbits  computed 

with  their  method  are  used  in  this  effort  to  define  the  nominal  path 
near  which  the  spacecraft  will  be  maintained. 

Trajectory  determination  for  a  spacecraft  that  moves  under  the 
influence  of  a.  two-body  system  of  forces  will,  however,  be  affected  by 
many  sources  of  error,  including  tracking  errors,  modeling  uncertainty, 
and,  possibly,  control  input  errors.  Orbit  determination  error 
analysis  seeks  to  quantify  the  impact  of  the  numerous  errors  that 

affect  the  motion  of  the  spacecraft.  The  result  of  the  error  analysis 
is  a  determination  of  the  spacecraft  position  and  velocity  uncertainty 
after  some  predetermined  period  of  flight  during  which  the  spacecraft 
is  affected  by  both  the  uncertainties  in  the  forces  and  the  errors  in 
tracking  data.  The  combined  magnitude  of  the  errors  may  be  found  to 
vary  depending  on  the  size  and  shape  of  the  spacecraft  orbit.  A 

reduction  in,  or  a  more  accurate  estimation  of,  the  magnitudes  of  the 

individual  errors  may  be  possible  and  could  then  lead  to  a  significant 

reduction  in  overall  vehicle  position  and  velocity  uncertainty. 

This  reduced  level  of  position  and  velocity  uncertainty  may,  in 
turn,  reduce  orbital  "maintenance"  costs,  such  as  the  propellant 
required  to  keep  the  spacecraft  near  the  nominal  orbit.  The  orbital 
maintenance  routine  is  referred  to  here  as  "station-keeping."  This 
cost  is,  in  part,  dependent  on  the  accuracy  of  the  tracking  information 
because  position  updates  using  inaccurate  tracking  data  may  result  in 
inefficient  use  of  control  energy  and  may  also  lead  to  unacceptable 
spacecraft  drift  from  the  nominal  path.  Other  error  sources  may  also 
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affect  spacecraft  drift  from  the  reference  trajectory  and,  therefore, 
may  increase  station-keeping  costs. 

This  research  is  concerned  with  quantifying  the  resulting  impact 
of  some  of  the  possible  error  sources  on  orbit  determination.  These 
investigations  use  a  nominal  path  for  the  spacecraft  orbit  that  is 
numerically  computed  as  a  solution  in  the  elliptic  restricted 
three-body  problem.  Simulated  tracking  data,  including  assumed  levels 
of  the  associated  errors,  can  then  be  produced.  The  overall  spacecraft 
position  and  velocity  uncertainties  can  be  computed  at  a  specified 
epoch  and  can  be  compared  as  functions  of  various  levels  of  error 
magnitudes  as  well  as  differing  orbital  shapes  and  sizes.  Using  two 
nominal  orbit  typos  and  various  input  parameter  errors  as  a  basis,  the 
resulting  errors  can  then  be  compared.  Chapter  1  briefly  covers  the 
background  of  the  elliptic  restricted  three-body  problem.  Chapter  2 
then  summarizes  throe  methods  of  orbit  determination  error  analysis, 
and,  finally.  Chapter  3  covers  the  results  obtained.  Follow-on 
research,  not  included  in  this  effort,  can  then  use  the  resulting  error 
levels  in  simulations  of  various  station-keeping  algorithms. 


CHAPTER  1:  BACKGROUND 


In  this  chapter,  the  elliptic  restricted  three-body  problem  and 
the  associated  coordinate  systems  are  reviewed;  the  equations  of  motion 
for  an  infinitesimal  mass  moving  in  the  gravity  fields  of  two  massive 
bodies  are  then  presented.  Next,  locations  of  the  libration  points  are 
discussed.  The  state  transition  matrix  and  the  construction  of  bounded 
nominal  orbits  near  the  collinear  Lagrange  points  are  then  summarized. 
Finally,  curve  fitting  the  nominal  trajectory  is  covered.  A  more 
thorough  discussion  of  these  topics  is  presented  in  Gordon. 


A.  Elliptic  Restricted  Three-Body  Problem 

The  elliptic  restricted  three-body  problem  is  a  simplification  of 
the  general  problem  of  three  bodies.  In  the  general  three-body 
problem,  each  of  the  three  bodies  is  assumed  to  be  a  particle  of  finite 
mass  and,  thus,  exerts  an  influence  on  the  motion  of  each  of  the  other 
bodies.  Neither  the  general  nor  the  restricted  problem  of  three  bodies 
has  a  general  closed-form  solution.  However,  when  problem 
simplifications  are  made,  particular  solutions  can  be  determined.  If 
the  mass  of  one  of  the  bodies  is  restricted  to  be  infinitesimal,  such 
that  it  does  not  affect  the  motion  of  the  other  two  massive  bodies 
(primaries),  the  restricted  three-body  model  results.  The  primaries 
are  assumed  to  be  in  known  elliptic  (ER3BP)  or  circular  (CT13BP)  orbits 
about  their  common  mass  center  (barycenter) .  The  problem  can  then  be 
completely  described  by  a  single  second-order  vector  differential 
equation  with  variables  appropriately  defined  for  a  specified 
coordinate  frame. 
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B.  Coordinate  Systems 


The  two  standard  coordinate  systems  used  in  the  analysis  of  this 
problem  have  a  common  origin  at  the  center  of  mass  (barycenter)  of  the 
primaries.  Primaries  with  masses  m^  and  m^  such  that  m^s  are 
assumed  here,  although  this  distinction  is  arbitrary.  The 

infinitesimal  mass  is  denoted  as  m  .  These  masses  (m  ,m  ,m  ) 

3  12  3 

correspond  to  particles  situated  at  points  P^,  P^,  and  P^, 

respectively.  The  barycenter  is  denoted  as  "B,"  and  the  resulting 

arrangement  is  shown  in  Figure  1-1.  The  rotating  coordinate  system  is 

defined  as  x  y  z  ,  and  the  , inertial  system  is  identified  as  x  y  z  . 

Note  that  both  coordinate  systems  are  right-handed,  and  the  x  and  y 

axes  for  both  systems  are  in  the  plane  of  motion  of  the  primaries.  The 

Xj  axis  is,  of  course,  assumed  to  be  oriented  in  some  fixed  direction; 

in  this  specific  formulation  of  the  problem,  it  is  assumed  to  be 

parallel  to  a  vector  defined  with  a  base  point  at  the  Sun  and  directed 

toward  periapsis  of  the  Earth’s  orbit.  The  rotating  x  axis  is  defined 

R 

along  the  line  that  joins  the  primaries  and  is  directed  from  the  larger 

toward  the  smaller  primary.  The  z  axes  are  coincident  and  are  directed 

parallel  to  the  primary  system  angular  momentum  vector.  The  y  axis 

R 

completes  the  right-handed  x  y  z  system. 

R  R  R 
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C.  Eiiuati'ons  of  Motion 


4 


Newtonian  mechanics  are  used  to  formulate  the  equations  of  motion 

for  in^  (the  spacecraft)  relative  to  B  as  observed  in  the  inertial 

reference  frame.  The  sum  of  the  forces  on  m^  resulting  from  both  the 

gravity  fields  of  masses  m  (the  Sun)  and  m  (the  Earth-Moon 

1  2 

barycehter)  and  from  the  solar  radiation  pressure  can  be  used  to 
produce  the  following  second-order  vector  differential  equation: 


„  m  m 

p  =  -  G  (— ^)  d  -  G  ( - ^)  r  + 

.3  3 

d  r 


(1-1) 


The  overbar  denotes  a  vector,  and  primes  indicate  differentiation 
with  respect  to  dimensional  time.  All  quantities  are  dimensional,  as 
appropriate,  and  the  quantity  "G"  is  the  universal  gravitiational 
constant.  The  scalars  "d"  and  "r"  in  equation  (1-1)  denote  the 
magnitudes  of  the  vectors  d  and  r,  respectively,  as  depicted  in 
Figure  1-1.  The  dimensionless  scalar  "k"  is  the  solar  reflectivity 
constant,  and  "S"  is  the  solar  radiation  pressure  constant.  The 
formulation  of  the  solar  radiation  force  model  and  the  values  for  the 
solar  radiation  constants  are  derived  from  previous  work  by  Bell. 

The  values  of  the  constants  are  described  in  Gordon. 

The  position  vector  p  is  written  in  rotating  components  as 


p  =  X 


X  +  y  y 


+  z  z 


(1-2) 


*  where  x^,y^,z^  are  unit  vectors.  The  velocity  and  the  acceleration  of 

the  spacecraft  (particle  with  mass  m^)  relative  to  the  barycenter  13 
«  as  observed  in  the  inertial  reference  frame  can  then  be  described.  The 

following  kinematic  expression  for  p"  can  be  derived: 
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p"  =  (x"-0"y-20'y'-0'^x)x  +(y"+0"x+2O'x'-0'^y)y  +z"z  .  (1-3) 

R  R  R 


Three  scaled  equations  of  motion  for  can  be  derived  using  the 
the  following  scaling  factors: 


(1)  The  sum  of  the  masses  of  the  primaries  equals  one 
mass  unit, 

(2)  The  mean  distance  between  the  primaries  equals  one 
unit  of  distance. 

(3)  The  universal  gravitational  constant  is  equal  to  one 
unit  by  proper  choice  of  characteristic  time. 

The  dimensional  equations  of  motion  can  be  simplified  and  scaled 
by  introducing  the  characteristic  quantities  defined  above  and  by 
introducing  the  nondimensional  mass  ratio  p,  "psuedo-potential"  U,  and 
the  scaled  solar  radiation  constant  s: 


and 


M  = 


(1-4) 


U  = 


(1-p) 

d 


+ 


0^(x^  +  y^ 


k  s 

■  "T 


(1-5) 


where  the  dot  denotes  the  derivative  with  respect  to  characteristic 
time.  The  scaled  solar  radiation  constant,  s,  is  derived  by  using  the 
characteristic  quantities  described  above.  Then,  the  vector  magnitudes, 
"d"  and  "r,"  are  written  in  terms  of  scaled  quantities  as: 
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•  ^  f  2^X/2 

d  =  [(x  +  M  R)  +  y  +  z  ]  , 


(1-6) 


r  =  ((x  -  R  +  fi  R)^  + 


(1-7) 


If  the  primaries  are  assumed  to  be  moving  in  a  circular  orbit, 
equations  (1-8),  (1-9),  and  (1-10)  reduce  to  three  scalar  equations  in 
the  simplified  form: 
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r» 


^  -  "ST  “  ^x’ 


(1-11) 


V  +  2  X  =  ~  =  U 

y  +  ^  ^  ay  ^y’ 


(1-12) 


au 

"ST 


=  u  . 


(1-13) 


The  scalar  equations  can  be  used  to  locate  the  five  libration 
points  in  the  rotating  reference  frame. 


D.  Locations  of  the  Lagrangian  Points 

By  using  scalar  equations  (1-11),  (1-12),  and  (1-13)  for  motion  in 
the  CR3BP,  the  locations  of  the  stationary  equilibrium  points  can  be 
determined.  Equations  (1-8),  (1-9),  and  (1-10)  can  be  used  to 
determine  ratios  of  distances  that  are  constant  in  the  ER3BP‘,  these 
ratios  are  related  to  the  locations  of  libration  points  that  have  been 
defined  in  the  ER3BP  and  that  "pulsate"  with  respect  to  the  rotating 
reference  frame  as  the  distance  between  the  primaries  varies  with  time. 


1.  the  CR3BP 

In  the  CR3BP,  the  five  libration  points  are  equilibrium  points  and 
are  stationary  with  respect  to  the  rotating  coordinate  frame,  that  is, 
they  are  locations  at  which  the  forces  on  the  third  body  sum  to  zero. 
The  arrangement  of  points  and  the  corresponding  nondimensional 
distances  are  depicted  in  Figure  1-2.  Note  that  three  of  the  libration 
points  (L  ,  L  ,  L  )  are  collinear  with  the  primaries;  one  collinear 

12  3 
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point  (L.  )  is  interior  to  the  primaries.  The  remaining  two  points  (L^ 
and  L^)  are  located  at  the  vertices  of  two  equilateral  triangles  that 
are  in  the  plane  of  primary  rotation  and  that  have  a  common  base 
between  the  primaries. 

In  the  CR3BP,  the  libration  points  are  stationary  in  the  rotating 
coordinate  frame.  Stationary  points  are  defined  as  points  at  which  the 
relative  velocity  and  acceleration  are  zero,  such  that 

x  =  y*z«x  =  y»z*0.  (1-14) 


By  using  equations  (1-14)  in  equations  (1-11)  through  (1-13),  the 
useful  conditions  U  =  U  =  U  =  0  are  found.  The  three  collinear 

X  y  z 

libration  points  can  be  readily  located  by  further  noting  that 
y  =  z  =0  for  the  points  located  on  the  rotating  x  axis. 

n 


2.  The  ER3BP 

Five  libration  points  also  exist  in  the  ER3BP,  but  they  are  not 

stationary  relative  to  the  rotating  frame;  rather,  the  collinear  points 

pulsate  along  the  x  axis,  and  the  triangular  points  pulsate  relative 

R 

to  both  the  x  and  the  y  axes  as  the  distance  between  the  primaries 
R  R 

varies  with  time.  The  equilibrium  solutions  can  be  located  by  using 
equations  (1-8)  through  (1-10)  to  find  ratios  of  certain  distances 
that  are,  in  fact,  constant  in  the  problem.  The  collinear  libration 
points  in  the  ER3BP  can  be  found  by  assuming  x  ^  0,  x  *  0,  and  y  =  y  = 
z=z=y=z=0.  The  relative  locations  of  the  libration  points  in 
the  ER3BP  are  depicted  in  Figure  1-3. 
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E.  State  Transition  Matrix 

The  state  transition  matrix  is  used  in  the  calculation  of  the 
acceptable  nominal  trajectory,  and  it  must  also  be  available  at  varying 
time  intervals  along  the  nominal  path  for  orbit  determination  error 
analysis  investigations  and  station-keeping  studies.  The  transition 
matrix  is  derived  in  connection  with  a  linearizing  analysis. 

The  equations  of  motion  for  the  infinitesimal  mass  in  the  ER3BP 
can  be  linearized  about  a  reference  trajectory  (nominal  path)  that  is  a 
solution  of  the  differential  equations.  The  states,  three  position  and 
three  velocity,  and  the  state  vector  x  are  defined  as 


X  = 
1 


X.  X,  »  y, 


X,  X, 


X  ® 
6 


Z, 


(1-15) 


and 


X  =  ^3. 


(1-16) 


The  reference  trajectory  is  defined  as  x^^^.  Therefore,  using  a 
Taylor’s  series  approach,  the  expansion  about  the  reference  path  is 
written  in  the  form  of  the  first-order  variational  equation 


d 

dt 


(  x)  =  X  =  A(t) 


(1-17) 


where  x  =  x  -  x  is  understood  to  be  the  vector  of  residuals  relative 
REF 

to  the  nominal  solution,  and  the  matrix  A(t)  contains  the  first-order 
terms  in  the  Taylor’s  series  expansion  of  the  equations  of  motion  about 
the  nominal  or  reference  solution  of  interest. 


1* 


Using  equations  (1-8)  through  (1-10),  A(t)  can  be  expressed  as 


A(t)  = 


0 


[Urr+GH 


I 

2en 


(1-18) 


where  all  four  submatrices  are  dimension  3x3  and 


with 


Uxx 

Uxy 

Ux* 

Urr  * 

Uyx 

Uyy 

Uyz 

Uxx 

Uzy 

Uxz 

(1-19) 


r  0 


1  o' 

0  0  . 
0  0 


In  equation  (1-19),  the  notation  is  simplified  for  the  partial 
derivatives;  for  instance 


3^U 


ax 


2 


Uxx. 


The  matrix  A(t)  can  then  be  evaluated  along  the  reference  trajectory. 

The  vector  differential  equation  (1-17)  governing  the  state 
variations  from  the  nominal  path  has  a  solution  of  the  form 
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(1-20) 


x(t)  =  #(1,1^)  x(t^) 


where  *(t,t  )  is  the  state  transition;  matrix  at  time  "t"  relative  to 
time  "t^."  The  matrix  then,  represents  the  sensitivities  of  the 
states  at  time  "t"  to  small  changes  in  the  initial  conditions.  It  is 
determined  by  numerically  integrating  the  matrix  differential  equation 


—  i(t,t  )  =  iCt.t  )  =  A(t)  tCt.t  ),  (1-21) 

dt  0  0  0 

with  initial  conditions  ♦(t  ,t  )  =  I,  the  6x6  identity  matrix.  Thus, 

0  0 

the  nonlinear  equations  of  motion  in  (1-8)  through  (1-10)  and  the 
matrix  equation  (1-21)  combine  to  result  in  42  first-order  differential 
equations  that  can  be  simultaneously  integrated  numerically  to 
determine  the  state  vector  and  its  associated  transition  matrix  at  any 
instant  of  time.  The  reference  trajectories  that  are  of  interest  in 
this  research  are  generated  by  a  numerical  integration  method  that  uses 
a  differential  corrections  process  developed  by  Howell  and 
Pernicka. The  orbits  include  solar  radiation  pressure  forces  as 
formulated  by  Bell^^^^  specifically  for  an  orbit  associated  with  the 
interior  Lagrange  point  in  the  Sun-Earth  system.  The  numerical 
integration  routines  used  in  this  work  are  fourth-  and  fifth-order 

[17] 

Runge-Kutta  formulas  available  in  the  386-Matlab  software  package. 


F.  Bounded  Orbits  Near  Libration  Points 

The  computation  of  bounded  periodic  and  quasi-periodic,  orbits  in 
the  vicinity  of  libration  points  has  been  of  increasing  interest  during 
the  past  100  years.  This  section  first  discusses  the  stability  of  the 
libration  points  in  the  CR3BP  and  the  ER3BP.  The  construction  of 
bounded  orbits  near  the  collinear  Lagrange  points  is  then  summarized. 
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Finally,  the  specific  reference  trajectories  used  in  the  orbit 
determination  error  analysis  and  station-keeping  studies  in  this  work 
are  introduced. 


1.  Stability  of  the  Libration  Points  in  the  CR3BP 

The  accomplishments  of  those  researchers  who  have  constructed 

bounded  orbits  near  collinear  libration  points  are  particularly 

significant  because  the  collinear  points  are  considered  "unstable" 

points  of  equilibrium  but  with  (only)  one  mode  producing  positive 

exponential  growth.  Bounded  motion  in  their  vicinity,  therefore,  is 

determined  by  deliberately  not  exciting  the  unstable  mode.  A  second 

mode  produces  negative  exponential  orbital  decay  and  is  also 

deliberately  not  excited.  In  the  CR3BP,  the  remaining  four  eigenvalues 

are  purely  imaginary.  The  existence  of  initial  conditions  that  result 

in  only  trigonometric  (sinusoidal)  functions  as  solutions  means  that 

the  collinear  libration  points,  while  unstable,  possess  conditional 

stability  (with  proper  choice  of  initial  conditions)  in  the  linear 
(181 

sense. 

The  triangular  libration  points  are  marginally  stable  in  the 
linear  sense  for  a  specific  range  of  primary  mass  ratio  in  the  CR3BP. 
Purely  imaginary  roots  in  two  conjugate  pairs  are  obtained  for  p:s.0385, 
which  is  given  here  to  four  decimal  places  and  is  sometimes  referred  to 
as  Routh's  value.  The  mass  ratios  (listed  here  to  three  decimal 

places),  for  example,  in  the  three-body  systems  of  the  Earth-Moon 
(|i  =  1.216  X  lO"*),  Sun-Ear th+Moon  (p  =  3.022  x  10~^)  and  Sun-Jupiter 

”4 

(p  =  9.485  X  10  )  all  satisfy  the  mass  ratio  requirement  for  marginal 
stability  of  the  triangular  points  in  the  linearized  model.  Natural 
satellites,  such  as  the  Trojan  asteroids  or  a  moon  of  Saturn,  occupy 
linearly  stable  orbits  near  triangular  libration  points  in  their 
respective  systems. 


2.  Stability  of  the  Libration  Points  in  the  ER3BP 


Several  researchers  have  analyzed  the  stability  of  the  libration 

points  in  the  elliptic  problem,  where  both  the  mass  ratio,  p,  and  the 

[18*'22] 

primary  orbit  eccentricity,  e,  influence  stability.  The 

instability  of  the  collinear  libration  points  as  determined  in  the 

circular  problem  for  all  the  values  of  mass  parameter  persists  for  the 

elliptic  problem:  an  analysis  of  the  collinear  points  shows  instability 

for  any  combination  of  the  values  of  both  p  and  e. 

The  results  of  a  linearized  stability  analysis  regarding  the 

effects  of  eccentricity  and  mass  ratio  on  the  linear  stability  of  the 

triangular  points  have  been  published  by  Danby  and  then  later  by 
[22] 

Bennett  .  Both  Danby  and  Bennett  have  numerically  generated  graphic 
depictions  of  the  linear  stability  region  in  the  p-e  plane.  For  the 
eccentricity  in  the  Sun-Ear th+Moon  ER3BP,  the  value  of  p  which  ensures 
linear  stability  is  only  slightly  less  than  Routh’s  value  (decreased  by 
approximately  one  percent).  An  interesting  aspect  of  the  p-e  stability 
region  is  that  a  range  of  values  of  p  greater  than  Routh’s  value  also 
defines  a  region  of  linear  stability  for  a  specific  range  of  values  of 
e  less  than  .3143. 

3,  Construction  of  Bounded  Collinear  Libration  Point  Orbits 

The  initial  goal  in  the  process  of  generating  bounded  orbits  near 
a  collinear  (unstable)  libration  point  is  to  avoid  exciting  the 
unstable  mode  associated  with  the  linearized  motion.  The  meteoric  dust 
particles  that  may  be  orbiting  near  Lagrange  point  in  the  Sun-Earth 
system  could  only  linger  near  that  point  if  they  arrive  with  the 
"correct"  initial  position  and  velocity  states  relative  to  L^.  The 
"correct"  initial  conditions  will  only  (primarily)  excite  the  stable 
modes  associated  with  the  linearized  motion  and  not  (or  minimally) 
excite  the  unstable  mode.  The  degree  to  which  the  unstable  mode  is 
excited  will  determine  the  length  of  time  that  the  dust  particles 
linger  near  L^. 
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The  third-order  analytic  representation  is  used  in  this  work  to 
provide  the  initial  model  for  the  trajectories.  The  method  of 
successive  approximations,  using  the  linearized  solution  as  the  first 
approximation  to  the  nonlinear  orbital  path,  and  the  method  of  dual 
time  scales  are  used  to  derive  the  third-order  result.  ’  ’  The 
method  of  successive  approximations  is  used  to  generate  an  asymptotic 
series  in  an  appropriately  small  parameter.  (The  square  root  of  the 
eccentricity  of  the  primary  orbit,  that  is  the  orbit  of  the  Earth-Moon 
barycenter  about  the  Sun,  is  the  small  parameter  used  here. )  The 
method  of  dual  time  scales  is  used  to  convert  the  system  of  ordinary 
differential  equations  into  a  system  of  partial  differential  equations. 
In  general,  the  method  of  multiple  scales  allows  the  various  nonlinear 
resonance  phenomena  to  be  included  in  the  approximate  analytic  solution 
and  provides  a  method  to  remove  secular  terms.  (Here,  "secular"  refers 
to  terms  that  include  the  time  variable  and  is  derived  from  the  French 
"si^cle"  meaning  century. ) 

[71 

The  analytic  solution  of  Richardson  and  Cary  for  the 
Sun-Ear th+Moon  ER3BP  has  been  derived  to  fourth  order,  but  the  third- 

[9*14 ] 

order  approximation  is  found  to  be  sufficient  for  this  research. 

A  numerical  integration  algorithm,  using  a  differential  corrections 
procedure  that  is  designed  to  adjust  the  first  guess  as  obtained  from 
the  analytic  approximation,  can  then  be  used  to  numerically  generate 
the  orbit  of  interest.  A  method  developed  by  Howell  and  Pernicka^®''^’^ 
is  used  here  to  generate  the  orbital  paths.  Their  method  initially 
employs  the  approximate  analytic  solution  to  compute  target  points.  A 
two-level  (position  matching  then  velocity  matching),  multi-step 
differential  corrections  algorithm  is  used  to  construct  a  numerically 
integrated,  bounded  trajectory  that  is  continuous  in  position  and 
velocity.  A  solar  radiation  pressure  model  developed  by  Bell*^^^  is 
also  incorporated  in  the  numerical  integration  procedure. 

The  method  of  Howell  and  Pernicka,  including  solar  radiation 
pressure,  uses  an  initial  analytic  guess  that  represents  a  halo  orbit 
or,  alternatively,  a  considerably  smaller  Lissajous  path.  The 
higher-order  terms  tend  to  slightly  alter  the  first-order  periodic  or 
quasi-periodic  path.  Consequently,  the  initial  target  path  for  a  halo 
orbit  will  generally  not  be  precisely  periodic.  The  two-level, 
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multi-step  differential  corrections  procedure  then  adjusts  the  initial 
analytic  target  orbit  and,  therefore,  will  compute  a  halo-type  orbit 
that  is  nearly  (but  not  exactly)  periodic. 


4>  The  Reference  Paths  Generated  for  This  Work 

Precisely  periodic  halo  orbits  exist  in  the  OtSBP.  They  also 
exist  in  the  ER3BP,  but,  in  the  ER3BP,  they  are  multiple  revolution 
trajectories  with  periods  much  longer  than  those  of  interest  here. 
Nearly  periodic  orbits  are  more  practical  in  the  ER3BP  and  are  much 
more  likely  to  be  used  in  mission  planning:  therefore,  the  goal  here 
should  be  slightly  modified  to  be  the  comparison  of  Lissajous  and 
"halo-type"  orbits.  The  general  shapes  of  the  three-dimensional 
halo-type  and  Lissajous  orbits  can  be  seen  by  plotting  three 
orthographic  views  of  each  orbit,  using  the  tabular  data  from  the 
numerical  integration  routine.  Figure  1-4  depicts  three  orthographic 
views  of  point  plots  for  the  Lissajous  orbit  used  in  this  research. 
Figure  1-5  contains  three  orthographic  views  (on  a  slightly  different 
scale)  of  the  considerably  larger  halo-type  orbit.  (Note  that,  in 
general,  the  amplitude  ratio  for  Lissajous  trajectories  is  arbitrary. 
In  halo  orbits,  however,  constraining  the  amplitude  ratio  results  in 
equalized  frequencies  for  in-plane  and  out-of-plane  motion. )  The 
orbits  are  depicted  in  the  rotating  reference  frame  centered  at  L^. 

Both  orbits  are  clearly  not  periodic;  a  Lissajous  orbit  is  often 
called  a  quasi -per iodic  path,  and  these  two  orbits  could  clearly  be 
termed  quasi-periodic  or  Lissajous  paths.  The  major  difference  between 
the  orbits  is  the  larger  size  of  the  halo-type  orbit;  however,  other 
differences  are  also  present.  The  maximum  x  and  y  excursions  of  the 
halo-type  orbit  are  approximately  four  times  as  large  as  those  of  the 
Lissajous  path.  Furthermore,  the  direction  of  motion  (clockwise  versus 
counterclockwise),  as  viewed  in  the  y-z  orthographic  depiction,  is 
different  for  the  two  orbits  used  here.  The  direction  of  motion  on  the 
halo-type  orbit  is  counterclockwise  in  the  y-z  depiction;  the  direction 
of  motion  is  clockwise  in  the  y-z  depiction  for  the  Lissajous  path. 
(Both  orbits  include  clockwise  motion  in  the  x-y  depiction.) 
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Figure  1-4.  Three  Orthographic  Views  of  a  Lissajous  Orbit. 
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Figure  1-5.  Three  Orthographic  Views  of  a  Halo-Type  Orbit. 
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The  two  orbits  can  also  be  differentiated  in  terms  of  Lhi.' 
direction  of  the  maximum  z  excursion  in  the  x-z  depiction.  If  the 
maximum  z  excursion  is  in  the  positive  z  direction,  the  orbit  can  be 
termed  a  member  of  a  "northern  family"  of  orbits.  When  the  maximum  z 
excursion  of  the  orbit  is  in  the  negative  z  direction,  the  orbit  is 
termed  a  member  of  a  "southern  family"  of  orbits.  In  the  x-z 
orthographic  depiction,  the  smaller  (Lissajous)  path  can  be  seen  to  be 
a  member  of  a  northern  family  of  orbits,  while  the  halo-type  orbit  is  a 
member  of  a  southern  family  of  orbits. 

Future  work  with  these  two  orbits  will  include  studies  that 
generally  require  access  to  a  nominal  path  that  is  at  least  piecewise 
smooth.  Some  method  of  curve  fitting  the  numerically  integrated  data 
must  consequently  be  investigated. 


5.  Curve  Fitting  the  Nominal  Path 

The  primary  goal  of  this  work  is  the  completion  of  orbit 
determination  error  analysis  for  libration  point  trajectories.  The 
conventional  method  for  solution  of  state  estimation  problems,  and  the 
technique  used  in  this  effort,  involves  linearizing  the  nonlinear 
equations  of  motion  about  a  reference  solution  (nominal  path)  and  then 
applying  linear  estimation  techniques.  The  orbit  determination  process 
is  thus  changed  from  estimating  the  state  of  a  nonlinear  system  to 
estimating  the  linear,  time-varying  deviations  from  the  reference 
trajectoi’y. 

The  reference  solution  used  in  this  research  is  generated  by 

numerical  integration  of  the  nonlinear  equations  of  motion.  In  one 

study,  an  investigation  that  used  a  consistent  dynamic  model  for  all 

(8) 

comparisons,  Richardson  has  shown  that  a  slight  reduction  in  fuel 
expenditure  can  be  realized  if  a  numerically  integrated,  rather  than  an 
approximate  analytic,  nominal  path  is  computed.  The  numerical 

[9-14] 

integration  method  developed  by  Howell  and  Pernicka  is  used  here 

to  generate  a  set  of  reference  points  for  both  position  (three  states) 
and  velocity  (three  states),  relative  to  the  libration  point  of 
interest,  at  specified  times.  Time  series  point  plots  of  all  six  state 
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variables  for  approximately  a  2-year  segment  of  a  Lissajous  orbit  are 
depicted  in  Figures  1-6  (position  states)  and  1-7  (velocity  states). 
The  method  computes  numerical  data  for  the  six  states  in  a  reference 
frame  that  is  centered  at  the  libration  point  (in  this  case  L^)  and 
that  rotates  with  the  primaries.  However,  the  state  estimation 
techniques  and  station-keeping  algorithms  used  in  this  work  require 
access  to  a  continuous  nominal  path,  rather  than  point  plots,  of 
acceptable  accuracy. 

The  reference  trajectory,  represented  as  a  (piecewise)  smooth 

curve,  could  be  constructed,  approximately  or  exactly,  through  the 

points  obtained  from  the  numerical  integration  routine.  The  work  here 

assumes  that  a  curve  that  passes  through  the  numerical  data  (exactly) 

is  preferred.  The  effort  required  to  generate  a  numerical  solution, 

including  forces  modeled  consistent  with  the  ER3BP  (or  even  more 

accurately  modeled  with  ephemeris  data)  would  seem  to  be  wasted  if  the 

reference  curve  deviates  too  far  from  the  numerical  data.  However,  a 

method  that  approximates  a  smooth  curve  through  the  points  is  also 

desirable;  that  is,  linear  interpolation  between  the  numerical  data 

[9] 

points  was  not  considered  acceptable.  In  one  study,  Pernicka  found 
that  station-keeping  costs  for  a  libration  point  orbit  were,  in  fact, 
sensitive  to  the  accuracy  of  the  curve  fit.  Clearly,  a  piecewise 
linear  curve  fit  could  not  accurately  match  the  concavity  of  the  actual 
orbital  path  between  data  points,  regardless  of  the  size  of  the  time 
steps  used  in  the  numerical  integration  routine. 

Four  methods  of  generating  a  curve  for  the  nominal  trajectory  have 

been  evaluated:  Fourier  series,  least  squares,  weighted  least  squares, 

and  cubic  splines.  The  states  associated  with  a  quasl-per iodic  path 

were  thought  to  be  the  most  difficult  to  curve  fit;  therefore,  various 

Lissajous  trajectories  were  used  to  evaluate  the  curve-fitting  methods. 

After  several  curve  fitting  evaluations,  a  cubic  spline  interpolation 

routine  was  selected  to  be  used  to  model  the  reference  trajectory  in 

the  state  estimation  simulations.  The  comparisons  of  the  curve  fit 

[15] 

methods  are  fully  described  in  Gordon. 
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CHAPTER  2:  ORBIT  DETERMINATION  ERROR  ANALYSIS 


Complete,  exact  knowledge  of  the  state  of  a  spacecraft  in  orbit  is 
generally  not  possible.  Individual  state  variables  cannot  be  measured 
precisely,  and  available  measurements  are  usually  some  function  of  the 
state  variables.  For  instance,  a  spacecraft  moving  along  a  libration 
point  trajectory  in  the  Sun-Earth  system  may  be  tracked  using  range  and 
range-rate  measurements  containing  random  errors.  The  spacecraft  may 
be  affected  by  forces  not  included  (or  inadequately  represented)  in  the 
dynamic  model,  and  model  parameters  may  be  uncertain.  By  definition, 
the  linearized  system  of  equations  used  to  model  the  nonlinear  state 
variations  is  a  further  approximation.  Also,  actual  control  inputs  may 
vary  slightly  in  magnitude  and  direction  from  those  commanded.  These 
sources  of  error  make  knowledge  of  the  spacecraft  state  uncertain. 
Computation  of  the  most  likely  current  state  of  the  spacecraft  in  the 
presence  of  measurement  and  model  uncertainty  is  the  focus  of  orbit 
determination. 

Error  analysis  involves  an  investigation  of  the  impact  of  various 
sources  of  error  on  orbit  determination.  The  output  of  an  error 
analysis,  as  used  in  this  work,  provides  the  magnitudes  of  state  vector 
variances  and  covariances,  thus  quantifying  the  relative  contributions 
of  the  significant  error  sources.  This  output  could  then  be  used  to 
predict  how  an  improvement  in  measurement  accuracy,  for  instance,  would 
lessen  state  uncertainty.  One  benefit  of  more  accurate  knowledge  of 
the  state  might  be  a  reduction  in  station-keeping  costs.  A 
mathematical  procedure  can  be  developed  to  combine  all  information 
about  the  spacecraft  state  and  filter  this  observed  data  based  on  the 
varying  degrees  of  uncertainty.  The  filter  then  produces  a  "best 
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estimate"  of  the  state  and  additionally  quantifies  the  resulting  state 
variable  uncertainties  in  preparation  for  an  error  analysis,. 

The  first  goal  of  this  research  is  to  investigate  orbit 
determination  error  analysis  for  libration  point  trajectories.  The 
conventional  method  for  solution  of  state  estimation  problems,  and  the 
technique  used  in  this  effort,  involves  linearizing  the  nonlinear 
equations  of  motion  about  a  reference  solution  (nominal  path)  and  then 
applying  linear  estimation  techniques.  The  orbit  determination  process 
is  thus  changed  from  estimating  the  state  of  a  nonlinear  system  to 
estimating  the  linear,  time-varying  deviations  from  the  reference 
trajectory.  The  reference  solution  used  here  is  generated  by  numerical 
integration  of  the  nonlinear  equations  of  motion.  The  integration 
routine  produces  a  tabular  listing  of  each  of  the  six  state  variables, 
indexed  by  the  time  steps  used  in  the  numerical  computation.  The 
reference  path  is  then  modeled  by  passing  a  cubic  spline,  with  the 
independent  variable  being  time.  Individually  through  each  of  the  six 
state  variables  (three  position  and  three  velocity) , 

This  chapter  includes  a  brief  background  summary  of  filtering 
methods,  and  the  three  methods  used  for  error  analysis  in  this  work  are 
then  discussed.  These  three  methods  are  the  Kalman  filter,  batch 
weighted  least  squares  covariance  analysis,  and  consider  covariance 
analysis.  Each  of  these  methods  computes  a  covariance  matrix  at  a 
specified  epoch,  and  the  diagonal  entries  in  the  covariance  matrix  are 
the  variances  of  the  state  vector  elements.  The  square  roots  of  these 
six  entries  are  estimates  of  the  state  uncertainty  levels  at  the 
indicated  epoch. 

The  Kalman  filter  is  fully  derived  in  the  discussion  below  for  use 
in  both  orbit  determination  and  error  analysis  investigations.  The 
remaining  two  methods  (covariance  analysis  and  consider  covariance 
analysis)  are  derived  only  for  error  analysis,  but  either  method  can  be 
used  as  an  important  step  in  the  derivation  of  a  batch  weighted  least 
squares  filter.  The  following  chapter  includes  a  discussion  and 
comparison  of  the  results  derived  from  implementation  of  these  three 
methods.  The  state  vector  variance  levels  determined  from  the  error 
analysis  studies  for  both  Lissajous  and  halo  orbits  can  then  be  used  in 
future  station-keeping  simulations. 
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A.  BACKGROUND 


The  historical  and  mathematical  background  associated  with 
filtering  theory  and  error  analysis  originates  with  least  squares 
theory.  The  least  squares  techniques  based  on  the  theories  of  Gauss, 
Legendre,  and  Adrain  were  powerful  methods  that  initially  were  used  to 
explain  the  motions  of  the  heavenly  bodies.  In  particular,  the 
development  of  a  sequential  algorithm,  recursive  least  squares,  by 

[24] 

Gauss  is  fundamental  to  modern  linear  filtering  theory.  Gauss  also 
selected  a  weighting  scheme  for  errors  that  allowed  his  weighted  least 
squares  approach  to  result  in  a  linear,^  unbiased,  minimum  variance 
estimate.  In  Great  Britain  during  1912,  Sir  Ronald  Aylmer 

Fisher  introduced  the  maximum  likelihood  method. The 
analysis  of  variance,  as  used  in  error  analysis,  was  also  developed 
chiefly  by  Fisher,  and  he  introduced  the  terms  "variance"  and  "analysis 
of  variance. 

In  the  years  surrounding  World  War  II,  the  driving  forces  behind 

the  continuing  research  involving  estimation  algorithms  changed.  The 

original  focus  was  motivated  by  the  need  to  track  objects  at  low 

altitudes  near  the  Earth.  During  World  War  II,  probabilistic  methods 

were  recognized  to  be  valuable  in  addressing  problems  in  radio 

communications,  radar,  and  anti-aircraft  defense.  In  the  early 

1940s,  Norbert  Wiener,  in  the  United  States,  and  Andrei  Nikolaevich 

Kolmogorov,  in  the  Soviet  Union,  developed  input-output  relationships 

for  statistically  optimal  f liters. Both  investigated  the 

problem  in  the  frequency  domain.  Wiener  provided  optimal  estimates 

only  in  the  steady-state  regime  for  a  continuous  system,  while 

Kolmogorov  also  treated  the  discrete-time  problem.  In  1942,  the 

first  explicit  solutions  for  least-squares  estimates  associated  with  a 

stochastic  process  were  provided  by  Wiener  using  stationary  signal  and 
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noise  processes.  Wiener  found  that  the  solution  of  the  Wiener-Hopf 
integral  equation  could  be  used  to  determine  the  covariance  function  of 
the  estimation  process.  He  approached  the  problem  by  using  Fourier 
analysis  and  developed  an  integral  equation  for  determining  the  impulse 
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response  matrix  for  the  filter.  However,  Wiener-Kolmogorov  filters 
and  related  fpllow-on  research  results  had  some  of  the  following 
drawbacks: 

(1)  They  were  complicated,  often  requiring  the  solution  of 
auxiliary  differential  and  algebraic  equations  and 
determination  of  the  roots  of  polynomials. 

(2)  They  Were  not  easily  updated  with  increases  in  the 
observation  interval. 

(3)  They  could  not  be  conveniently  adapted  to  a  vector 
approach. 

After  World  War  II,  much  of  the  research  was  redirected  to  the  problem 
of  orbital  motion  (but  now  of  spacecraft,  not  planets).  With  the  new 
focus,  implementation  required  that  the  filtering  methods  be  adaptable 
to  vectorized  models. 

In  the  1950s  and  after,  research  into  determination  of  satellite 
orbits  was  a  driving  force  in  the  development  of  state  space  filters. 
The  Apollo  program  established  the  need  and  provided  the  resources  for 
development  and  application  of  optimal  filters  and  also  focused  the 
efforts  of  many  in  the  scientific  community. The  frequency  domain 
integral  methods  of  Wiener  and  Kolmogorov  gave  way  to  the  methods  of 
least  squares  adapted  to  dynamic  models  that  were  described  by 

differential  equations.  J.W.  Follin  at  Johns  Hopkins  University  in 
1955  and  M.  Blum  in  1958  developed  recursive  approaches  based  on  this 
idea  and  laid  the  foundation  for  the  development  of  dynamic 
filters. In  1959,  Peter  Swerling  published  a  method  of 

127  301 

estimating  satellite  orbits  near  Earth  from  observational  data.  ’ 

He  developed  an  innovations  process  for  updating  position  estimates. 

Swerling’ s  work  was  the  first  known  attempt  to  tackle  updates  in  a 

vector  estimation  problem,  and  his  useful  recursive  formulation  can  be 

specialized  to  obtain  the  Kalman  filter. 

Rudolf  Emil  Kalman  is  usually  credited  with  the  development  of  a 

state-space,  nonstationary,  recursive  filter  that  is  in  wide  use  today. 

However,  the  claim  is  not  universally  accepted.  Swerling 

unsuccessfully  challenged  the  editors  of  the  AIAA  Journal  to  give  him 

[271 

credit  for  the  discovery  of  the  Kalman-type  filter.  Swerling’ s 
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update  equations  for  the  error  covariance  were  in  a  much  more 
cumbersome  form  than  those  produced  by  Kalman;  and  Kalman’s  algorithm, 
though  more  restrictive  than  Swerling’s,  was  particularly  suited  to  the 
dynamic  state  estimation  problems  associated  with  spaceflight 
research.  In  March  1960,  Kalman  published  his  first  description  of 

[341 

the  recursive  filter.  In  April  of  1960,  Richard  H.  Battin 

independently  developed  a  recursive  optimal  filter  similar  to  that  of 

roc  ocl 

Kalman.  ’  Then,  in  November  1960,  R.L.  Stratonovich  independently 

obtained  the  Kalman  filter  equations  by  treating  the  linear  filter  as  a 

special  case  in  nonlinear  filtering. 

Initially,  Kalman  derived  and  published  the  discrete  model  of  his 

filter.  In  1961,  Kalman  and  Richard  S.  Bucy  published  a  continuous 

version  of  the  sequential  filter,  commonly  called  the  Kalman-Bucy 
[30] 

filter.  Research  interest  in  the  Kalman  filter  continued  with 

Stanley  F.  Schmidt’s  development  of  a  Kalman  filter  that  could  be 
implemented  using  an  on-board  digital  computer  for  guidance  in  a 
circumlunar  mission.  The  linear  theory  that  forms  the  basis  of  the 
Kalman  filter  could  have  been  a  handicap  for  such  mission 
Implementation;  however,  Schmidt  linearized  the  equations  of  motion 
about  the  nominal  path  and  formulated  the  filter  to  update  (extend)  the 
nominal  path  after  every  position  fix.  Schmidt  first  applied  this 
extended  filter  to  the  C-5A  aircraft  navigation  system. 
Consequently,  the  extended  Kalman  filter  is  sometimes  identified  as  the 
Kalraan-Schmidt  filter. 


B.  The  Kalman  Filter 

[38  39] 

Kalman’s  early  published  articles  ’  were  somewhat  difficult 
to  understand.  When  Schmidt*^®*  initially  tried  to  apply  an 
operational  Kalman  filter,  he  and  his  engineering  staff  found  it 
necessary  to  have  initial  and  follow-up  briefings  from  Kalman.  Schmidt 
found  the  notation  and  the  concepts  used  in  the  filter’s  development  Lo 
be  difficult  to  grasp.  Strang*^^^  recently  suggested  that  the 

sections  of  most  texts  that  include  discussions  of  the  Kalman  filter 
are  rather  discouraging  to  read.  Every  author  has  his  or  her  own 
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notation;  even  with  similar  notation,  one  Kalman  filter  can  look 
different  from  another  because  authors  also  have  several  choices 
between*  equivalent  equations  for  propagation  of  the  error  covariance 
ands  Kalman  gain  matrices.  So,  even  beyond  notational  differences, 
propagation  equations^  may  appear  to  vary  in  a  mathematical  development 
of  the  Kalman  filter. 

The  Kalman  filter  is  a  least  squares  formulation  that  incorporates 
[261 

a  dynamic  model.  For  a  spacecraft  in  orbit,  the  governing  dynamic 
model,  that  is,  the  set  of  equations  of  motion,  is.  combined  with 
available  measurements  to  improve  knowledge  of  the  spacecraft  state. 
The  linearized  model  and  the  state  transition  matrix  used  in  the  orbit 
determination  investigations  have  been  previously  described.  The 
measurement  process  is  added  in  this  section,  and  the  equations  are 
written  for  the  measurement  residuals  relative  to  the  nominal  path. 
Recall  that  the  nominal  path  is  defined  by  six  cubic  splines  passing 
through  the  numerical  data  representing  the  states.  Then,  the 
measurement  and  dynam*'  models  (and  the  associated  notation)  used  in 
the  following  derivation  are 


Measurements:  z  =M  x  +  v  (2-1) 

k  k  k  k 


Dynamic  Model:  x  =  ♦(k+l,k)  x  ,  (2-2) 

k^-l  k 


where  z  is  the  measurement  residual  vector  at  time  step  k;  x  is  the 

k  k 

residual  state  vector  at  time  step  k;  M^  is  the  measurement  matrix  that 

is  linearized  about  the  nominal  path;  ♦(k+l,k)  is  the  state  transition 

matrix  at  time  step  k+1  relative  to  time  step  k;  and  v^  is  the 

measurement  noise  vector  and  is  assumed  to  have  the  statistics 
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(2-4) 


E(v  v^)= 

k  k 


V 

k 


where  "E"  is  the  expectation  operator,  0  is  the  zero  vector <  and  is 
the,  measurement  noise  covariance  matrix.  (Process  noise  is 
incorporated  in  this  model  in  a  later  section;  control  inputs  are  not 
considered  in  these  orbit  determination  error  analysis  studies. ) 

In  this  work,  both  range  and  range-rate  tracking  are  simulated. 
The  measurement  matrix,  M^,  is  linearized  about  the  nominal  path 
represented  by  six  cubic  splines  passing  through  the  numerical  data  for 
the  six  spacecraft  states.  The  two  scalar  nonlinear  tracking  equations 
are: 


range  =  R  =  ((x-x  )^+  (y-y  +  (z-z 

■  s  a 


(2-5) 


range-rate  =  R  =  [(x-x  )(x-x  )+(y-y  ) (y-y  )+(z-z  )(z-z  ))/R.  (2-6) 

a  a  a  a  a  a 


Here,  the  subscript  "s"  denotes  the  coordinates  associated  with  the 
tracking  station.  The  measurement  matrix,  M,  is  a  time-varying  matrix, 
derived  from  equations  (2-5)  and  (2-6)  through  a  linearizing  process, 
and  is  written 


M  = 


m  m  m  0  0  0 

11  12  13 


m  m  m  ro  m  m 

21  22  23  24  25  26 


(2-7) 


where 
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(2-8) 


m  =  m  =  (x-x  )/R  , 

11  24  8 


"12“  "25'  (y-y.)'Ti  . 


(2-9) 


m  =  m  =  (z-z  )/R  , 

13  26  a 


(2-10) 


m  =  [(x-x  )R^-NUM(x-x  )]/R^. 

21  •  a 


(2-11) 


"'22'  ((y-y,)R^-NUM(y-ya)]/R^, 


(2-12) 


m  =  Kz-z  )R^-NUM(z-z  )1/R^, 

23  8  a 


(2-13) 


NUM  =  [(x-x  )(x-x  )  +  (y-y  )(y-y  )  +  (z-z  )(z-z  )]. 

a  a  8  8  a  a 


(2-14) 


If  M  represents  the  matrix  M  evaluated  at  time  t  ,  the  time-varying, 

Ic  ^ 

linearized  vector  measurement  equation  at  time  t^  is  then  given  by 


=  M  X 


k  k 


+  V 

k 


(2-15) 
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for  residuals  from  the  nominal  path.  The  time-varying  measurement 

matrix,  M^,  has  the  nonzero  entries  (2-8)  through  (2-13)  that  are 

evaluated  along  the  nominal  path.  The  measurement  residual  vector,  z^, 
represents  measured  deviations  from  the  nominal  path.  The  vector  v^, 
as  the  measurement  noise  vector,  indicates  that  the  updates  have  some 
level  of  uncertainty. 

A  measurement  that  arrives  at  step  k  and  the  system  propagation 
equation  that  provides  an  estimate  of  the  state  at  step  k  can  be 

combined  using  weighted  least  squares.  (See  Gordon^'^^  for  a  complete 

discussion  of  the  least  squares  and  weighted  least  squares 

approaches. )  The  system  equation  can  be  propagated  to  step  k  from  the 
"best"  estimate  of  position  at  step  k-1  using  the  state  transition 
matrix.  One  key  to  determining  how  to  weigh  (filter)  these  two 

estimates  (the  measurement  and  the  state  propagated  from  previously 
filtered  measurements)  is  to  use  the  covariance  matrix  in  a  weighted 
least  squares  process.  The  covariance  matrix  that  depends  on  the 
measurements  through  step  k-1  can  be  denoted  by  Pj^  j-  The  covariance 
matrix  at  step  k,  before  the  measurement  at  step  k  is  incorporated,  is 
denoted  by  P..  Propagation  of  the  state  covariance  matrix  can  be 

illustrated  by  first  defining  x  as  the  estimate  propagated  to  state  k 

.V  ^  ...  .V 

before  measurement  z  ,  x  as  the  estimate  after  measurement  z  ,  and  x 

k  k  k  k 

as  the  true  state  propagated  to  state  k. 

The  propagation  equations  can  then  be  written 


X  =  *(k,k-l)  X  , 
k  k-1 


(2-16) 


X  =  ♦(k,k-l )  X 
k  k-1 


(2-17) 


P  =  E[(x 

k-1  k-1 


A  ^ 

-X  ,)(x 

k-1  k- 


A  .T, 
-X  )  J  , 
1  k-1 


(2-18) 
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P  =  E[(x  -  X  )(x  -  X  )’^1, 
k  k  k  k  k 

=  E[(*(k,k-1)^  -♦(k.k-l)x^  ^M^Ck.k-Dx^  -*(k,k-l)x^  j”^] 

■>  k-l  k-1  k-1  k-1 

=  E[*(k,k-l)(x  -X  )(x  -X  )V(k,k-l)] 

k-1  k-1  k-1  k-1 

=  ♦(k,k-l)  E[(^  -x^  )(x^  -x^  )'^]*'^(k,k-l).  (2-19) 

k-1  k-1  k-1  k-1 

Hence,  using  (2-18)  in  (2-19),  an  equation  for  mapping  the  covariance 
matrix  from  state  k-1  to  state  k  is  given  by: 

P ,  =  ♦(k,k-l)  P  t'^(k,k-l)..  (2-20) 

k  k-1 

The  measurement  residual  vector  z  can  be  combined  with  the  state 

k 

estimate  at  step  k  in  a  weighted  least  squares  solution  by  utilizing 

additional  definitions.  For  instance,  the  state  uncertainty  vector  at 

time  t  is  defined  as  c  and 
k  k 

X  ’=  X  +  c  (2-21) 

k  k  k 

where  the  state  uncertainty  is  assumed  to  have  the  statistics 

E(c  )  =  0  (the  zero  vector),  (2-22) 

k 
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(2-23) 
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The  matrices  f2-25),  (2-26),  and  {2-21)  can  be  used  to  combine 

equations  (2-1)  and  (2-21)  into  a  matrix  equation 


A  X  +  e  =  b. 

k 


(2-29) 


The  weighted  least  squares  formulation  for  this  matrix  equation  uses 
the  inverse  of  the  matrix  V,  defined  in  (2-28),  for  weighting  the 
errors  in  an  unbiased,  minimum  variance  scheme.  Minimizing  the  cost 
function  J(x),  where  J(x)  *  e\~^e,  is  completed  by  substituting  (2-29) 
into  the  cost  function  and  then  setting  the  first  derivative  (of  the 
cost  function)  with  respect  to  x  equal  to  zero.  This  computation 
returns  the  solution  for  the  optimal  state  estimate  x^,  which  minimizes 

the  cost  function  J(x),  at  time  t  as 

k 


Ox  (aV‘a)”^  aV*F. 

k 


(2-30) 


The  equations  (2-25)  through  (2-28)  can  then  be  used  in  equation  (2-30) 

i.  j  1411 
to  produce 


*M 


4.  P 


.  K-l: 


(2-31) 


This  solution  for  in  (2-31)  requires  the  inversion  of  a  matrix  that 
may  be  quite  large.  Kalman  developed  a  sequential  algorithm  that 
requires  inversion  of  a  smaller  matrix  at  each  step.  From  the 
definition  of  the  covariance  matrix  in  (2-18)  and  by  then  using  (2-30) 
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arid  (2-29),  in  that  order,  the  covariance  matrix  can  be  represented  as 


P  =  (aV*A)"V  (2-32) 

k 

Using  (2-26)  and  (2-28)  in  (2-32),  a  computationally  efficient  equation 

A 

for  P  can  be  derived: 

P  =  (mV^M  +  p■^)“^  (2-33) 

k  kkkk 

A 

While  the  covariance  i.  ix  P^  is  the  essential  matrix  used  to  complete 
the  error  analysis  (cq‘  i  .-nee  analysis)  portion  of  this  work,  the  full 
Kalman  filter,  includi..ij  the  state  vector  update  equation,  is  derived 
here  for  completeness.  The  full  Kalman  filter  may  also  be  quite  useful 
as  an  optimal  observer  in  follow-on  research  concerning  this  problem. 

The  Schurr  identity  or  "inside-out  rule"  can  be  used  to  give  a 
smaller  dimension  matrix  for  inversion: 


P 


k 


=  (mVX  +  P"*)"^ 

k  k  k  k 


=  P  -  P  m”^[  M  P  m’’  +  V  ]"^  M  P  . 

k  kk  kkk  k  kk 


(2-34) 


The  Kalman  gain  matrix  can  be  derived  from  (2-34)  as 


PytM  P  +  V  ]"V 
kkkkk  k 


(2-35) 
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Using  equation  (2-35) j  equation  (2-34)  can  be  rewritten  as 


P 

k 


K  M  P 

k  k  k 


II  -  K  M  ] 

k  k 


(2-36) 


This  allows  equation  (2-31)  to  be  expressed  in  a  more  computationally 
efficient  form: 


A 

X. 


II  -  K  M  1  P  (m3  S  +  P"' 

k  k  k  k  k  k  k 


=  [I  -  K  M  1 

k  k 


X 


tl  -  KM] 

k  k 


P.  M 


-1  ~ 
k  k 


(2-37) 


This  expression  for  x^  can  be  reduced  further  by  showing  that 


(I-  K  M  ]  P  V"*  »  K  ,  (2-38) 

k  k  k  k  k  k 

SO  that 

X  =  X  +  K  iS  -  X J.  (2-39) 

k  k  k  k  k  k 

The  relationships  in  equations  (2-16),  (2-20),  (2-35),  (2-36),  and 
(2-39)  can  be  combined  to  produce  a  weighted  least  squares  estimate  of 
the  (dynamic)  state  of  the  spacecraft  relative  to  the  nominal  path  and 
in  the  presence  of  measurement  errors.  These  five  equations  are  used 
here  to  determine  the  spacecraft  residual  state  at  a  given  epoch.  The 
covariance  matrix,  P^,  reflects  the  uncertainty  associated  with  a  given 
residual  state  at  the  selected  time.  A  final  adjustment  to  the  error 
covariance  equations  can  also  ensure  that  the  filter  accepts  at  least  a 
small  portion  of  each  measurement  residual  update. 
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One  problem  associated  with  the  filter  is  divergence;  one  reason 
for  this  divergence  can  be  the  filter's  failure  to  accept  measurement 
updates.  As  in  recursive  least  squares,  the  filter  can  "go  to  sleep" 
and  not  accept  any  substantial  part  of  new  measurements.  Similarly,  in 
batch  least  squares,  a  single  observation  has  a  decreasing  effect  on 
the  solution  as  the  number  of  observations  increases.  The  Kalman  gain 
matrix,  K  ,  depends  in  part  on  the  state  covariance  matrix,  as  defined 

^  T  —1 

in  equation  (2-33).  The  term  MV  M  will  be,  at  least,  positive 

..  k  k  k 

semi-definite;  thus,  may  decrease  as  measurements  continue  to 
arrive,  and  the  gain  matrix  can  be  rewritten  using  equations  (2-36)  and 
(2-38)  as 


K  =P  m”^  V‘\  (2-/10) 

k  k  k  k 

A 

Clearly,  if  P  decreases,  then  K  will  decrease  as  well.  The 

k  ^  k 

measurement  updates,  Iz  ~M  x  1,  in  equation  (2-39)  will  be  given  little 

K  K  K  ^ 

weight,  and  the  update  equation  (2-39)  will,  in  effect,  become  x^^  =  x^^, 
ignoring  measurements.  Modeling  errors  will  cause  the  system  to  drift 
from  Xj^,  and  measurements  to  correct  the  state  will  not  be 
incorporated.  One  solution  to  the  divergence  problem  is  incorporation 
of  system  error  covariance  in  the  Kalman  f liter. This 
approach  is  reasonable  because  the  linear  model  is,  in  fact,  an 
approximation,  and  parameter  errors  and  neglected  forces  add  to  the 
system  uncertainty.  (However,  the  appropriate  magnitude  of  the  added 
system  uncertainty  is  somewhat  unclear  and,  perhaps,  arbitrary. ) 
Adding  system  noise  modifies  the  equations  as  follows: 


X  =  <Kk,k-l)  X  +  w 

k  k-l  k 


(2-41  ) 


where  w^  is  the  added  system  noise  vector  and  is  the  system  noise 
covariance  matrix  such  that 
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(2-42) 


E(w  )  =  0 

k 


E(w5^)  =  W^. 

k  k  k 


The  covariance  matrix  propagation  expression  (2-20)  then  appears  as 


P  =  *(k.k-l)  ,  ♦^(k.k-1)  +  W^. 

k  k-1  k 


(2-43) 


The  magnitudes  of  the  positive  entries  in  determine  the  relative 
weight  assigned  to  the  predicted  position  x^^  as  compared  to  the 
measurement  z^.  The  Kalman  filter  with  system  uncertainty  is  then 
defined  in  terms  of  equations  (2-16),  (2-43)  in  place  of  (2-20),  (2-35), 
(2-36),  and  (2-39). 

The  Kalman  filter  can  be  used  independently  to  perform  a 
preliminary  error  analysis.*^®*  The  computation  of  the  state 

A 

covariance  matrix  after  the  measurement  update  (P^)  provides  an 
indication  of  the  accuracy  of  the  state  estimate  at  a  specified  time 
step.  The  filter  equations  provide  a  relatively  straightforward  method 
of  generating  the  state  covariance  matrix.  When  a  nominal  trajectory 
is  available,  the  Kalman  gain  and  the  state  covariance  matrices  can  be 
precomputed  and  stored.  That  is,  only  the  measurement  covariance 
matrix  is  required  for  this  computation;  measurements  are  not  needed  to 
precompute  the  covariance  matrices. 

At  any  time,  the  diagonal  entries  in  the  state  covariance  matrix 
are  the  variances  of  the  individual  states,  and  the  off-diagonal 
elements  indicate  the  cross-correlation  between  the  states.  Also,  the 
square  root  of  the  sum  of  the  first  three  diagonal  entries  of  the 
covariance  matrix  (a  partial  trace)  is  an  indicator  of  the  overall 
uncertainty  in  the  position  states.  The  overall  velocity  state 
uncertainty  can  likewise  be  investigated.  In  general,  at  a  specified 
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epoch,  the  square  roots  of  the  diagonal  elements  of  the  covariance 
matrix  are  indicators  of  the  consequences  of  accumulated  uncertainties 
over  the  tracking  period. 

Several  methods  can  be  used  to  compare  error  contributions.  The 
error  budget  is  a  listing  of  individual  error  sources  and  the  relative 
contribution  of  these  sources  to  the  total  state  uncertainty.  A 
sensitivity  analysis  that  uses  various  levels  of  measurement  accuracy 
in  the  computations  can,  in  turn,  inform  the  filter  designer  of 
individual  error  source  contributions.  A  source-by-source  sensitivity 
analysis  may  be  valuable  in  assessing  potential  hardware  improvements 
or  measurement  scheduling  changes.  Sensitivity  curves  are  plots  of 
changes  in  the  uncertainty  of  an  individual  state  or  group  of  states  as 
measurement  accuracy  levels  are  varied.  Other  useful  graphical  output 
Includes  plots  of  state  uncertainty  (square  roots  of  the  diagonal 
entries  in  the  state  covariance  matrix)  and  plots  of  elements  of  the 
Kalman  gain  matrix  as  functions  of  time.  Both  types  of  plots  are 
depicted  later  in  the  results  section  of  this  chapter.  With  added 
system  noise,  these  plots  will  appear  to  approach  non-zero  steady-state 
values;  thus,  the  gains  do  not  approach  zero,  and  the  filter  does  not 
"go  to  sleep." 

C.  Batch  Weighted  Least  Squares 

Error  analysis  using  batch  weighted  least  squares  is  similar  in 
many  ways  to  that  of  the  Kalman  filter.  The  algorithm  for  computation 
of  the  covariance  matrix  in  batch  weighted  least  squares  is  derived 
here,  and  much  of  the  notation  and  derivation,  as  used  in  this  work,  is 
in  common  with  the  Kalman  filter.  Equations  (2-1)  through  (2-33)  from 
the  previously-described  Kalman  filter  formulation  are  also  used  (in  a 
similar  sequence)  in  the  derivation  of  batch  weighted  least  squares, 
and  these  equations  are  not  repeated  here. 

The  error  analysis  using  batch  weighted  least  squares  is  developed 
here  as  a  covariance  analysis.  (The  derivation  of  state  vector  updates 
is  straightforward,  but  is  not  necessary  in  this  work. )  For  this  error 
analysis,  the  diagonal  entries  of  the  covariance  matrix  after  a 
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predetermined,  but  variable,  number  of  measurement  updates  can  be  used 
to  predict  state  uncertainty  levels  used  in  the  station-keeping 
simulations  of  chapter  three. 

The  differences  between  the  Kalman  and  batch  weighted  least 
squares  filters  are  predominantly  in  the  combination,  or  filtering,  of 
position  updates.  In  the  Kalman  filter  formulation,  measurements  and 
the  propagated  state  estimate  are  combined  in  a  recursive  weighted 
least  squares  algorithm,  using  covariance  matrices  to  determine 
relative  weights.  The  updated  state  estimate  is  then  propagated 
forward,  and  the  process  of  combining  measurements  with  state  estimates 
is  continued  until  the  end  of  the  tracking  period  is  reached. 
Alternatively,  in  batch  weighted  least  squares,  numerous  measurement 
observations  are  treated  as  a  batch  of  position  updates  that  are 
reduced  to  a  single  epoch  and  combined  with  the  predicted  state  vector. 
Here,  that  epoch  is  denoted,  arbitrarily,  as  t^. 

In  order  to  facilitate  a  simple  example,  we  suppose  that  two 
measurements,  one  at  time  t^  and  one  at  time  t^,  need  to  be  combined 
with  the  dynamic  model’s  position  prediction  at  time  t^^.  Using 
equations  (2-25),  (2-26),  and  (2-29)  with  the  error  vector  e  omitted 
for  simplicity,  the  following  matrix  equation  can  be  formulated: 


•M 

b  = 


z 

1 

<v 

z 

2 

= 

= 

> 

X 

0 

X 

0 

I 

X  =  A  X  . 
0  0 


(2-44) 


In  general,  many  more  than  two  observations  will  be  combined  to  form 

T  —t  -1 

the  covariance  matrix  (A  V  A)  .  For  just  the  two  observations  in  the 
example  above, 
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(2-45) 


V  = 


For  this  example,  the  matrix  V  is  dimension  10x10.  The 

measurement  noise  matrices  V  and  V  are  each  dimension  2x2, 

1  2 

corresponding  to  uncertainties  in  the  range  and  range-rate 

measurements,  and  the  state  covariance  matrix  P  is  of  dimension  6x6. 

0 

A  reasonable  tracking  schedule  of  9  range  and  range-rate  measurements 

per  day  for  20  days  would  generate  a  matrix  V  of  dimension  366x366. 

T  -1  -1 

Inversion  of  the  matrix  V  in  the  covariance  computation  (A  V  A) 
could  then  lead  to  ill-conditioning  and  associated  computational 

difficulties.  In  one  study  for  ISEE-3,  the  error  analysis  simulations 

using  batch  least  squares  did  encounter  such  numerical 

difficulties. 

An  adjustment  to  the  weighted  batch  least  squares  algorithm 
that  permits  inversion  of  a  smaller  dimension  matrix  after  each 

measurement  update  could  reduce  possible  numerical  difficulties.  The 
matrix  elements  in  the  covariance  matrix  (a\”^A)"'  can  be  multiplied 
using  definitions  (2-44)  and  (2-45)  to  give 


P'^=  P"'+  *^(1  .t  )mV^M  ♦(!  ,t  )+  ♦'^(t  ,t  )mV^M  ♦(!  ,t  ).  (2-46) 

0  0  1  0  1  l  1  10  2’  0  2  2  2  2’  0 


The  computation  can  be  generalized  to  k  measurements  taken  after 

time  t^  to  give  an  equation  that  can  be  used  for  calculation  of  the 

state  covariance  matrix  P  at  epoch  t  : 

0  0 


=  P 


j:*^(tj.t^)M\ 

1=1 


M  $(t  ,t  ). 
1  110 


(2-47) 
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The  state  covariance  can  then  be  propagated  from  t^  to  t^^  by  using 
equatibh  (2-20)  to  give: 


K  =  (2-48) 

k  k  0  0  k  .0 

The  propagation  equation  (2-48)  incorporates  the  assumption  that 

A 

the  most  recent  measurement  update  arrived  at  t^^,  and  thus  is 
computed  rather  than  The  diagonal  entries  of  this  covariance 

matrix  are  the  variances  of  the  states  at  time  t^,  after  measurement 
updates  through  time  t^.  Notice  that  only  measurements  and  the  initial 
state  vector  (by  using  the  covariance  matrix  P^)  are  considered 
uncertain  in  this  formulation.  The  results  of  covariance  matrix 
computations  for  the  Kalman  filter  and  for  batch  weighted  least  squares 
will  be  compared  using  identical  tracking  schedules  and  statistics  in 
the  last  section  of  this  chapter.  However,  a  third  method  of 
covariance  matrix  calculation  that  incorporates  additional  sources  of 
error  is  first  derived. 


0.  Consider  Covariance  Analysis 

Consider  covariance  analysis  uses  a  batch  weighted  least  squares 
formulation  but  also  includes  parameter  uncertainty.  Model 
parameters  that  are  considered  uncertain  in  this  work  are  the  planetary 
masses  (through  the  dimensionless  mass  parameter  p),  the  locations  of 
the  tracking  stations,  and  the  solar  reflectivity  constant.  This 
derivation  has  some  formulas  in  common  with  the  Kalman  filter  and  with 
batch  weighted  least  squares  covariance  analysis;  however,  the  consider 
covariance  matrix  derivation  is  somewhat  more  complicated  than  both  of 
these  previously  described  methods.  So,  a  brief  overview,  including  a 
few  general  steps  in  the  problem,  is  provided  here,  before  the  full 
derivation  is  completed. 

The  goal  of  a  covariance  analysis,  as  used  in  this  work,  is  to 
compute  the  covariance  matrix  after  some  period  of  tracking.  The  state 
uncertainties  are  described  by  the  magnitudes  of  the  variances  found  on 
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the  diagonal  of  this  covariance  matrix.  In  general,  at  the  epoch  of 
interest.,  the  state  uncertainty  is  thought  to  be  the  consequence  of  the 
accumulated  uncertainties  in  the  model  and  the  measurements.  Some 
parameters  may  also  be  considered  uncertain. 

The  covariance  matrix  of  interest  is  assumed  to  be  the  aposteriori 

A 

(after  measurements  through  time  t  )  covariance  matrix  P  at  epoch  t  . 

k  k  k 

This  matrix  is  defined  in  equation  (2-18)  as 


=  E[(x  -X  )(x  -X 
k  k  k  k  k 


The  state  vector  denoted  by  x  is  the  true  (yet  unknown)  state;  the 
state  vector  denoted  by  x  is  the  "best"  estimate  of  the  true  state  at 

k 

time  t  .  Here  the  word  "best"  means  that  the  vector  x  might  be 
k  k 

computed,  for  instance,  by  minimizing  a  given  cost  function  appropriate 
to  the  problem  at  hand.  The  covariance  matrix  may  be  computed 
recursively,  as  in  the  Kalman  filter,  or  in  a  batch  formulation,  as  in 
batch  least  squares.  The  batch  formulation  is  used  here,  where  the 
measurements  are  reduced  to  an  initial  epoch  and  then  propagated  by 
using  the  transition  matrices.  The  addition  of  uncertain  parameters 
complicates  the  dynamic  and  measurement  models  and,  consequently,  the 

A 

computation  and  propagation  of  the  covariance  matrix  P  . 

[41 1  ** 

The  dynamic  model  is  now  written 


X  =  A  X  +  B  c. 


(2-49) 


where  the  matrices  A  and  B  are  computed  from  the  usual  linearizing 
process  relative  to  the  nominal  path  and  by  using  nominal  values  for 
all  parameters  considered  uncertain.  The  vector  c  represents  the 
residuals  of  the  parameters  considered  uncertain  relative  to 
predetermined  nominal  values.  Here, 


c  =  (k  p  X  y 

S  S 


(2-50) 
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where  "k"  is  the  dimensionless  solar  reflectivity  constant,  “/i"  is  the 
dimensionless  mass  parameter,  and  "x^,  y^,  z^"  are  the  scaled  tracking 
station  coordinates.  The  matrix  B  is  time  varying  and  is  given  by 


B  = 


0 
0 
0 
b 

41  42 

3  b 

SI  S2 

3  b 

61  62 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


(2-51) 


where 


s(x+Rfi) 


(2-52) 


sy 

b  =  — , 

SI  .3’ 
d 


(2-53) 


sz 

b  =  — , 

61  .3 

d 


(2-54) 


(x+Rp)-R(l-fi)+skR  pR+x-R(l-p) 


b  = 
42 


3(l-ji)  (x+Rp)  V3sk(x+pR)^  3Rp(x-R(l-p) )‘ 


(2-55) 


y  3(l-p) (x+Rp)yR-3sk(x+pR)yR  3Rfi(x-R(l-p) )y 

—  +  -  +  - , 

3  .5  s 

r  d  r 


(2-56) 


7. 


b  = 
62 


d 


3 


z  3(l-p)  (x+Rji)zR-3sk(x+fiR)zR  3R/i(x-R(l-fi)  )z 


(2-57) 
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The  statistics  associated  with  both  x  (the  true  slate  vector)  and 
c  (the  vector  whose  elements  are  the  true  parameter  values)  can  now  be 
described.  The  error  vector  (sometimes  called  the  state  uncertainty 
vector)  at  time  t^  is  defined,  as  in  the  Kalman  filter  and  batch  least 
squares  derivations,  as  and 


X  =  X  +  c 
k  k  k 


(2-58) 


where  the  state  uncertainty  is  assumed  to  have  the  statistics 


E(c  )  =  0 

k 


(2-59) 


E(c  ?)  =  P  . 

k  k  k 


(2-60) 


The  consider  parameter  error  vector  at  time  t^  is  defined  as  and 


c  =  c  +  8  , 
k  k  '  k 


(2-61) 


with  the  consider  parameter  uncertainty  assumed  to  have  the  statistics 


E(8  )  =  0. 
k 


(2-62) 


E(8  c^)  =  0  (the  zero  matrix), 
k  k 


(2-63) 


E(5  (33)=P^  =  P  • 

k  k  k  c 


(2-64) 
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The  matrix  is  the  consider  parameter  covariance  matrix  and  is 
assumed  in  this  work  to  be  stationary. 

The  observation  equation,  for  time  t^  is  now  given  as 

z  +m‘'c  +  v.  (2-65) 

k  k  k  k  k  k 

t 

The  matrix  is  equal  to  the  matrix  described  in  the  Kalman  filter 
derivation;  for  the  computation  of  the  consider  parameters  are 

assumed  to  be  at  their  nominal  values.  The  matrix  is  derived,  also 

k 

by  linearizing  the  nonlinear  measurement  equations,  while  assuming  that 
the  consider  parameters  vary  about  their  nominal  values  and  that  the 


state  variables  are  equal 

to 

their  nominal 

values. 

Both  matrices 

are  time  varying.  This  matrix 

is  given  by 

0 

0 

m  m  m 

11  12  13 

=  - 
If 

» 

(2-66) 

K 

0 

0 

m  m  m 

21  22  23 

where  the  six  nonzero  elements  of  are  equal  in  magnitude  but 
opposite  in  sign  to  the  indicated  elements  of  M*.  For  example,  the  1,3 
entry  in  the  matrix  is  equal,  as  indicated  in  equation  (2-66),  to 
The  statistics  associated  with  the  measurement  equation  are 


E(v  )  =  0, 

k 

E(c  v^)  =  0, 

k  k 


E(6  v^)  =  0, 

k  k 


E(v  v^)  =  V  =  V. 

k  k  k 


(2-67) 

(2-68) 

(2-69) 

(2-70) 
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In  the  remainixig  steps  of  this  derivation,  the  initial  time  is 
assumed  to  be  t^.  (Recall  that,  in  the  batch  weighted  least  squares 
derivation,  t^  was  chosen  as  the  initial  epoch — the  choice  is 
completely  arbitrary  in  the  derivation  and  t^^  was  chosen  here  to 
broaden  the  presentation. )  Therefore,  observations,  in  addition  to 
can  be  made  at  t  and  t  ,  for  example,  and  these  take  the  form: 

k+l  k+2 


z  =M  X  +M  c  +v, 
k>l  k4^1  kn  k-fl  M  k+l 


(2-71) 


2  =M  X  +M  c  +v. 
k*2  k*2  k*2  k*2  k*2  k*2 


(2-72) 


As  in  batch  weighted  least  squares,  these  observations  (through 
t  ,  for  instance)  can  be  reduced  to  a  specific  epoch,  such  as  t  . 

ktj  k 

However,  with  consider  covariance  analysis,  the  uncertain  parameters 

also  have  a  transition  matrix  0,  and  a  state  propagation  equation  such 

as  equation  (2-2)  cannot  be  used.  The  state  propagation  now  includes 

the  effects  of  variations  in  the  consider  parameters,  and  the  state 

[41] 

propagation  equation  is  given  by 


x(t  )  =  $(t  ,t  )x  +  9(t  ,t  )c  , 
k+l  k+l  k  k  k+l  k  k 


(2-73) 


i(t  .t  )  =  A(t)  $(t  , 

k+l  k  k+l 


(2-74) 


with 


and 


0(t  .t  )  =  A(t)  e(t  .t  )  +  B(t) 

k-*-!  k  k-fl  k 


(2-75) 


with 


e(t  .t  )  =  0. 

k  k 

where  0  denotes  the  zero  matrix  (here  it  is  6x5).  In  this  work,  we 

assume  that  the  parameter  uncertainty  does  not  vary  with  time; 

therefore,  c  =  c.  The  observations  for  the  k,  k+1,  k+2,  etc  time 
k 

steps  can  then  be  combined  to  give 

S  =  Mx  +  mS  +  7,  (2-76) 

X  c 

where 


M 

X 


M 


K  ) 

k't’l  k'fl  k 

M’'  ♦(t  ,t  ) 
k+2  k+2  k 


(2-77) 


M 


M" 


M**  0(t  , 

k+1  k+1 

M’'  0(t  , 

k+2  k+2 


t  )  +  M 

k 

t  )  +  M 

k 


c 

k+1 


c 

k+2 


(2-78) 
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The  consider  covariance  matrix  P  ,  after  a  predetermined  number  of 

C 

measurement  updates,  can.  be  computed  by  using  the  equation  (as  in 
•Weighted  or  unweighted  least  squares) 


P 

C 


=  p®  =  E[  (x  -X  )  (x  -X  )^] . 
k  k  k  k  k 


(2-79) 


Before  deriving  a  solution  for  (2-79)  in  terms  of  measurement  and 
model  (including  parameter)  uncertainties,  a  brief  overview  of  the 
remaining  computational  process  is  included.  This  consider  covariance 

A 

matrix,  P*^,  is  the  covariance  matrix  that  contains  the  consequences  of 

k 

model  (including  parameter)  and  tracking  uncertainty  over  the  specified 
tracking  period,  for  example  from  t^  through  order  to  compute 

this  covariance,  the  measurements  over  the  tracking  period  of,  perhaps, 
20  days  are  reduced  to  the  initial  epoch  of  t  by  using  equations 

'‘c  * 

(2-77)  and  (2-78).  The  matrix  P^  will  then  be  computed  and  propagated 
to  the  end  of  the  tracking  period,  which  is  denoted  here  as  t  .  The 

A 

propagated  coyariance  matrix  P®^^  contains  the  state  variances  along 
the  diagonal.  The  positive  square  roots  of  these  variances  will  be  the 
state  vector  element  uncertainties  (standard  deviations),  as  computed 
in  this  error  analysis  study. 

The  optimal  residual  state  estimate  x^^  is  required  in  (2-79)  and 
is  computed  through  a  weighted  least  squares  approach;  however,  when 
some  parameters  are  considered  uncertain,  a  more  complicated  derivation 

A  ^ 

of  X  (and  thus  P  )  needs  to  be  addressed.  As  in  the  Kalman  filter  and 

^  ~  -1— 
batch  least  squares  derivations,  the  cost  function  J(x)  =  e  V  e  is 

minimized  to  find  the  optimal  state  estimate  x.  The  vector  e  is 

defined  as  in  (2-27);  however,  the  expanded  cost  function  will  now  also 

include  consider  parameter  terms: 


J(x)  =  (z  -  M  X  -  M  c)  V 


^(z  -M  X  -  M  c)+(x-x)P  ^(x-x). 


(2-80) 
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After  setting,  the  first  derivative  of  (2-80)  equal  to  zero  and  after 
some  algebraic  simplification,  a  formula  for  x  (an  optimal  estimate  for 
x)  is  obtained: 


X  =  (M  V  M  +  P  )  ‘(M  V  (z  -  M  c)  +  P  x). 

XX  X  k  c 


(2-81) 


Clearly,  a  batch  or  recursive  algorithm,  using  a  reformulation  of 
(2-81),  could  be  developed.  This  algorithm  could  then  be  used  to 
predict  the  "best"  estimate  of  the  spacecraft  state,  where  the 

sa  fa 

propagated  state  (x)  and  the  measurement  (z)  are  appropriately  weighted 
by  using  the  inverses  of  the  covariance  matrices  P  and  V,  respectively. 
However,  the  focus  of  this  research  is  orbit  determination  error 
analysis,  and  only  the  aposteriori  covariance  matrix  must  be  computed 
to  complete  this  analysis.  Optimal  state  estimates  need  not  be 
computed.  For  convenience,  we  simplify  the  notation  by  letting 


p  =  (mV‘m  +  P"‘)"‘. 

X  X 


(2-82) 


Notice  that  P  is  identical  to  the  covariance  matrix  derived  in  the 
batch  weighted  least  squares  filter.  The  adjusted  observation  equation 


z  =  Mx  +  Mc  +  v  (2-83) 

k  X  c 


includes  the  assumption  of  uncertainty  in  c.  This  assumption  is  valid 
because  the  consider  parameter  errors  are  assumed  to  be  zero-mean,  as 
depicted  in  equation  (2-59);  hence,  E(c)  =  E(c).  Equation  (2-83)  is 
then  substituted  into  (2-81),  and,  by  also  using  definition  (2-82), 
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equation  (2-81)  simplifies  to 


X  =  P(mV^M  )x 

X  k  X 


+  P(mV^)(MP  +  v)  +  PP’^x, 

X  k  c 


(2-84) 


where  equation  (2-61)  is  also  used  in  the  simplification.  In 

A 

preparation  for  computation  of  the  aposteriori  covariance  matrix  P  , 

A  ~  ^ 

the  term  (x-x)  can  be  found  in  a  compact  form  as 


X  -  X  =  PP“‘e  +  P(mV‘)(M  P  +  v). 

X  k  c 


(2-85) 


Equation  (2-85)  is  computed  by  using  a  form  of  (2-82)  given  by 


P“^  =  (p-^-  mV^m  ) , 

X  k  X 


(2-86) 


and  by  also  using  equation  (2-58).  After  using  equations  (2-85), 

(2-70),  and  (2-64)  in  definition  (2-79)  and  after  some  algebra,  a 

^  rs6i 

compact  equation  for  P  is  derived: 

C 


P 

c 


=  P  +  PM  V"  M  P  M  V"  M  P  =  P  . 

X  c  c  c  X  k 


(2-87) 


This  is  the  state  consider  covariance  matrix  at  time  t^,  the  epoch 
to  which  the  measurements,  compiled  within  the  tracking  period,  are 
reduced.  Clearly,  the  consider  covariance  matrix  in  equation  (2-87) 
computes  a  covariance  matrix  at  least  as  large  as  the  covariance  matrix 

A 

(P)  computed  in  batch  weighted  least  squares  because  the  added  term 
(PM\  P  mV  V  P)  is  at  least  positive  semi-definite.  As  in  batch 

X  C  C  C  X 

weighted  least  squares,  this  covariance  matrix  must  be  propagated  from 
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some  initial  epoch  to  a  time  corresponding  to  the  end  of  the  tracking 
period.  the  end  of  the  tracking  period  is  selected  here  to  be  t 

k+J 

The  existence  of  the  consider  transition  matrix  complicates  the 

[41  ] 

covariance  propagation  equation,  where 


X  -  X 
k+j  k+J 


=  tCt  ,  t  ) (x  - 
k+j’  k  k 


0(t  .,t  )(c  -  c), 
k+j  k 


(2-88) 


E[(x  -  X  )(x  i-  X  )^). 

k+j  k+j  k+j  k+j 


(2-89) 


If  the  uncertain  parameters  were  not  considered,  equation  (2-89) 
could  be  simply  computed  by  using  a  form  of  propagation  equation 
(2-19).  However,  with  parameter  uncertainty,  expression  (2-88)  is 
substituted  into  equation  (2-89),  and  the  following  propagation  formula 
for  the  state  covariance  (after  some  algebra)  results: 


K  .  =  .-tJ  P"  ..t  )  +  e(t  .t  )  P“  e‘^(t  .t  )  + 

k+j  k+j  k  k  k+j  k  k+j  k  k  k+j  k 


♦  (t  ,t  )  P  V"^  M  P°  0'^(t  ,t  )  + 

k+j  k  k  X  c  k  k+j  k 


,.t  )  P^  V"^  M  ?  ♦^(t  ..t  ).  (2-90) 

k+j  k  k  c  X  k  k+j  k 


The  covariance  matrix  P°  can  also  be  denoted  here  as  P° 

k+j  k+j 

because  measurements  are  assumed  to  be  taken  through  time  t  .  This 

k+j 

aposteriori  covariance  matrix,  propagated  to  time  t  ,  contains  the 
state  variances  along  its  diagonal.  These  variances  are  the 
consequences  of  the  measurement,  parameter,  and  force  uncertainties 
that  have  accumulated  over  the  tracking  period  from  t^  until  t^  The 
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positive  square  roots  of  these  variances  can  then  be  used  as  the  state 
uncertainties  (state  standard  deviations).  The  results  of  this  error 
analysis,  or  more  specifically,  of  this  consider  covariance  analysis, 
as  compared  with  the  previous  two  methods  used  in  this  work  for  error 
analysis,  will  be  summarized  in  the  next  chapter. 
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CHAPTER  3:  ORBIT  DETERMINATION  ERROR  ANALYSIS  RESULTS 


How  do  the  results  incorporating  the  error  analysis  methods 
covered  in  this  chapter  compare  to  results  found  in  other,  similar 
efforts?  This  question  will  be  addressed  shortly,  but  first  a  survey 
of  the  input  error  levels  used  in  similar  error  analysis  studies  serves 
as  a  valuable  introduction.  It  is  noted  that  these  initial  uncertainty 
levels — for  instance,  measurement,  parameter,  and  force  uncertainties-- 
differ,  sometimes  significantly,  in  magnitude,  but  they  have  been  used 
in  several  orbit  determination  error  analysis  investigations.  The 
values  of  these  uncertainties  may  be  used  to  compute  the  diagonal 

entries  in  the  input  covariance  matrices  (such  as  V,  P^,  and  P).  The 

investigations  included  in  this  first  section  deal  with  libration  point 
halo  or  halo-type  orbits  near  in  the  Sun-Ear th+Moon  three-body 
system.  (Smaller  Lissajous  orbits  are  discussed  in  the  comparisons 
included  in  the  later  sections  of  this  chapter. ) 

In  the  next  three  sections,  error  analysis  results  that  are 
produced  from  each  of  the  three  methods  (Kalman  filter,  batch  weighted 
least  squares,  and  consider  covariance  analysis)  discussed  previously 
are  also  presented;  additionally,  the  results  corresponding  to 

selection  of  both  halo-type  and  Lissajous  trajectories  as  nominal 
orbits  are  compared.  The  error  analysis  results  from  several  other 
investigations  are  then  compared  with  the  consider  covariance  results 
found  in  this  effort.  The  error  levels  to  be  used  in  the 

station-keeping  simulations  of  chapter  three  are  then  computed. 
Finally,  in  the  last  results  section,  a  Monte  Carlo  simulation  is  used 
to  complete  statistical  hypothesis  tests  concerning  the  probability 
distributions  and  the  popuration  means  of  the  residual  state  errors. 
These  error  levels,  with  assumed  probability  distributions  and  means, 
can  then  be  used  in  Monte  Carlo  simulations  of  station-keeping 
algorithms  developed  in  future  research. 
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A.  Survey  of  Input  Error  Levels  Used  in  Other  Research 


Table  1  lists  the  input  error  levels  assumed  in  other  error 
analysis  studies  as  compared  to  the  levels  used  in  this  study.  The 
errors  are  denoted  by  the  symbols  generally  used  in  the  derivation 
sections  of  this  chapter.  The  dimensionless  solar  reflectivity 
constant  is  k;  the  tracking  site  location  uncertainty  is  S  and  is  input 
as  an  equal  uncertainty  level  for  each  of  the  site  coordinates  x  ,  y  , 

8  S 

z^;  range  tracking  is  R;  range-rate  tracking  is  RR;  and  the  uncertain 
mass  parameters  are  p  for  Earth,  p  for  the  Sun,  and  p  for  the  Moon. 

6  8  m 

The  last  column  contains  the  uncertainty  in  the  dimensionless  mass 
parameter  p  that  would  be  "equivalent"  to  the  errors  listed  for  the 
individual  mass  parameters.  (Recall  that  p  =  (p  +  p  )/(p  +  p  +  p  )  for 

on  8  6  n 

the  three-body  system  of  interest  in  this  work. )  The  approximate  value 
of  <r^  (the  standard  deviation  of  p)  is  calculated  from  extensive 
(10,000)  Monte  Carlo  trials  while  assuming  the  individual  mass 
parameter  error  levels  that  are  described  in  each  of  these  studies.  It 
is  this  value  of  mass  parameter  uncertainty  (o*^)  that  is  used  in  the 
formulation  of  the  consider  covariance  analysis  used  in  this  work. 

An  entry  in  Table  1  that  is  " — "  means  that  the  study  in  question 
did  not  clearly  indicate  if  an  uncertainty  of  this  type  was  used.  The 
error  levels  listed  for  one  entry  (Longuski^^®^ )  are  computed  from 
currently  used  mission  planning  data  provided  during  private 
conversations  with  the  several  technical  experts  listed  for  this  source 
in  the  references  section.  Some  studies^^^"^**^  also  use  additional 
input  error  sources  in  their  orbit  determination  error  analysis 
simulations;  however,  the  errors  used  in  this  work  include  those  found 
to  be  most  significant  in  the  literature.  Consequently,  the  error 
levels  identified  in  Table  1  for  this  investigation  were  then  input  as 
the  standard  deviations  (squared,  as  appropriate,  for  the  variances 
along  the  diagonals  of  the  measurement  and  parameter  covariance 
matrices)  for  the  covariance  analysis  studies  summarized  in  the 
following  sections. 
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Table  1.  Survey  of  Error  Analysis  Input  Errors. 


ONE  STANDARD  DEVIATION  ERRORS 


STUDY 

k 

4 

s 

(km) 

R 

(km) 

RR 

(m/sec) 

M  M 

e  • 

(  <-  km^/sec^ 

^  ) 

cr 

M 

(x  p) 

Mistretta^^^’ 

15*/. 

.010 

.007 

1.000 

3.08 

xlO® 

.0726 

Joyce 

10*/. 

.002 

.015 

.002 

.3986 

1.327 

xlO^ 

.0490 

1.411 

xlO"”^ 

Efron'^^*^^' 

Rodriguez- 

10*/. 

.002 

.015 

.002 

3986 

1.327 

-xlO^ 

.0490 

1.411 

xlO”'^ 

Canabal ^ 

— 

.010 

.015 

.003 

— 

— 

“ 

Longuski 

13% 

.0003 

.010 

.001 

.4903 

4030.7 

.0100 

1.231 

xlO"’^ 

This  Work 

13% 

.010 

.015 

.003 

.3986 

1.327 

xlO* 

.0490 

1.411 

xlO"’^ 
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B;  folman  Filter  Results 


The  Kalman  filter  generates  an  aposteriori  state  covariance  matrix 
that  can  readily  be  used  for  error  analysis  studies.  It  also  computes 
optimal  gain  matrices  that  are  used  to  provide  a  relative  weight  to  the 
current  measurement  versus  a  prediction  of  the  actual  spacecraft 
position.  Therefore,  in  addition  to  covariance  analysis,  the  Kalman 
filter  may  prove  useful  as  an  optimal  tracking  filter  in  future 
research.  Its  tracking  abilities  are  only  partly  introduced  here  for  a 
spacecraft  following  a  nominal  halo-type  orbit.  The  results  of  error 
analysis  using  the  Kalman  filter  are  then  included  in  the  next  section, 
where  the  results  obtained  from  using  the  Kalman  filter  and  weighted 
batch  least  squares  filter  covariance  analyses  are  compared. 

Results  from  simulations  incorporating  the  Kalman  filter  described 
previously  show  both  good  tracking  and  the  anticipated  gain  propagation 
characteristics.  These  characteristics  can  be  depicted  in  several 
ways.  Figure  3-1  is  a  plot  for  the  residual  y(t)  state  estimate  and  is 
produced  from  the  results  of  one  Monte  Carlo  simulation  that  was 
continued  for  30  days.  Plots  for  the  other  five  state  variables  show 
similar  behavior,  but  they  are  not  included  here.  The  full  state 
measurement  residual  vector  is  corrupted  by  random  noise  and  is  input 
to  the  filter.  The  filter  output  should  "follow"  this  noisy 
measurement  signal  but,  at  the  same  time,  be  consistent  with  the 
propagation  of  the  state  residual.  The  positive  and  negative  values  of 
the  square  root  of  the  (2,2)  element  of  the  state  covariance  matrix, 
plotted  as  functions  of  time,  form  a  one  standard  deviation  envelope 
for  filter  uncertainty.  This  envelope  approaches  the  zero  error  level 
because  this  formulation  of  the  Kalman  filter  does  not  use  added  system 
noise.  These  solid  (standard  deviation)  lines  appear  in  Figure  3-1 
just  above  and  below  the  horizontal  level  that  corresponds  to  zero 

errors  in  the  residual  y(t)  state.  The  line  " . "  depicts  the 

tracking  signal,  including  errors  in  the  y(t)  residual  state,  that  is 
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The  line " . "  is  the  y  state  output  from  the  kalman  filter 

Figure  3-1 .  Kalman  Filter  Plots  for  the  y  Stale  Variable. 
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input  for  the  Monte  Carlo  run.  The  dashed  line  is  the  Kalman 

filter  output.  Nine  "noisy"  measurements  are  simulated  per  day,  and 
the  filter  tends  to  smooth  out  tracking  uncertainty  by  continuing  to 
decrease  the  weighting  that  the  measurements  are  given  in  the  recursive 
least  squares  solution. 

The  output  depicted  in  Figure  3-1  is  from  a  Kalman  filter 
formulated  using  equations  (2-16),  (2-20),  (2-35),  (2-36),  and  (2-39). 
When  system  uncertainty  is  not  added  to  the  formulation,  measurements 
are  given  less  weight  (compared  to  the  propagated  state)  as  the 
simulation  progresses.  The  lines  depicting  the  standard  deviation 
envelope  in  Figure  3-1  can  be  seen  to  continue  to  approach  the 
zero-error  level,  and,  as  this  standard  deviation  decreases,  the 
measurements  can  also  be  seen  to  have  a  decreasing  effect  on  the  filter 
output.  That  is,  the  filter  output  can  be  seen  to  "smooth  out"  and 
approach  the  zero-error  level  as  measurements  continue  to  arrive.  The 
filter  would  eventually  ignore  measurement  updates,  and  this  tendency 
of  the  Kalman  filter  is  often  called  "going  to  sleep."  This 

formulation  of  the  filter  (without  added  system  noise)  is  used  in  the 
covariance  analysis  simulations  described  in  this  section  because  the 
arbitrary  level  of  added  system  noise  would  otherwise  affect  the 
covariance  analysis:  however,  to  improve  the  tracking  capabilities  of 
the  filter,  added  system  noise  would  be  used  in  a  formulation  of  an 
optimal  tracking  filter  that  could  be  used  in  future  research. 

Figure  3-2  is  a  plot  of  the  residual  y(t)  state  estimate  and  is 
produced  from  the  results  of  one  Monte  Carlo  simulation  of  a  Kalman 
filter  that  includes  added  system  noise  in  its  formulation.  The  Kalman 
filter  incorporating  added  system  noise  uses  equations  (2-16),  (2-35), 
(2-36),  (2-39),  and  (2-43);  the  addition  of  the  system  noise  matrix 
in  equation  (2-43)  allows  both  the  predicted  spacecraft  position  and 
the  measurement  updates  to  be  given  continuing,  relative  weights  in  the 
sequential  weighted  least  squares  formulation  of  the  Kalman  filter. 
The  standard  deviation  envelope,  simulated  measurements,  and  the  filter 
output  are  depicted  in  Figure  3-2  with  the  same  characterizations  used 
in  Figure  3-1.  The  one  standard  deviation  envelope  in  this  formulation 
decreases  to  a  steady-state  level  within  1  day.  The  level  of  added 
system  uncertainty  is  slightly  exaggerated  in  this  example  in  order  to 
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The  line " . "  is  the  y  state  output  from  the  kalman  filter 


Figure  3-2.  Kalman  Filter  Plots  for  the  y  State  Variable  when  Added  System  Noise  is  Used. 
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more  clearly  depict  the  differences  evident  in  the  two  types  of  Kalman 

filter  derivations.  The  filter  output  for  the  y  state  ("-  - ")  tends 

to  "follow"  the  noisy  measurement  information  (" . ")  throughout  the 

30  days  of  tracking  more  closely  than  the  filter  output  presented  in 
Figure  3-1  tends  to  follow  the  "noisy"  measurements.  The  addition  of 
system  uncertainty,  in  effect,  decreases  the  relative  level  of 
confidence  given  numerically  to  the  propagated  state  in  the  filter  and 
provides  a  consistent  weight  in  the  least  squares  computation  for  the 
measurement  updates. 

Figures  3-3  and  3-4  show  the  Kalman  gain  matrix  element  (2,2)  for 
two  different  formulations  of  the  filter.  Figure  3-3  depicts  the 
decrease  in  the  gain  matrix  element  for  a  Kalman  filter  that  does  not 
include  added  system  noise.  The  magnitude  of  this  entry  decreases 
toward  a  zero  value  as  measurements  continue  to  arrive  in  the  30-day 
tracking  period.  Note  that  the  gain  matrix  K^^  in  equation  (2-39) 
determines  the  relative  weight  that  measurements  are  given  in  the 
filter  update  equation;  as  the  gain  elements  decrease,  the  filter  will 
tend  to  more  soundly  "go  to  sleep."  Figure  3-4  depicts  the  gain 
element  (2,2)  that  decreases  to  a  steady-state  value  after 
approximately  1  day  of  tracking.  (Recall  that  the  added  system  noise 
level  in  this  example  is  relatively  large  in  order  to  facilitate 
illustration  of  the  differences  between  these  two  formulations.)  The 
steady-state  gain  value  is  nonzero  because  this  Kalman  filter 
derivation  incorporates  added  system  noise;  a  nonzero  gain  value  allows 
measurements  to  continue  to  be  given  some  weight  in  the  filter’s  state 
propagation  equation.  Other  gain  matrix  elements  exhibit  similar 
characteristics,  but  they  are  not  shown  here. 

The  Kalman  filter  derivation  used  in  this  effort  for  orbit 
determination  error  analysis  (covariance  analysis)  incorporates  the 
error  sources  of  initial  state  uncertainty  (the  matrix  P^)  and 
measurement  noise  (the  matrix  V  ).  The  addition  of  system  uncertainty 

k 

(the  matrix  W^)  improves  filter  tracking  but  would  not  be  helpful  in 
completing  a  covariance  analysis.  Another  method,  batch  weighted  least 
squares  covariance  analysis,  is  similarly  derived  (to  a  point)  and  can 
be  used  to  possibly  confirm  the  results  obtained  using  the  Kalman 
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Gain  element  2,2 


Figure  3-3.  Gain  Element  2,2 
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Gain  element  2,2 


Figure  3-4.  Gain  Element  2,2  When  System  Noise  is  included  in  the  Filter. 
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filter.  In  the  following  section,  the  covariance  analysis  results 
using  a  formulation  of  the  Kalman  filter  that  does  not  include  added 
system  noise  (uncertainty)  is  compared  to  the  results  obtained  from  a 
weighted  batch  least  squares  filter. 


C.  Comparison  of  Weighted  Batch  Least  Squares  and  Kalman  Filters 

The  derivations  of  both  the  weighted  batch  least  squares  and  the 
Kalman  filters  are  somewhat  similar.  Assuming  use  of  the  same  level  of 
tracking  errors  and  initial  state  errors,  these  two  filters  should 
provide  the  same  stats  vector  uncertainty  levels  after  identical 
tracking  periods.  (The  agreement  may  be  approximate  due  to  numerical 

round-off  error. )  The  following  results,  in  Table  2,  correspond  to  180 
range  and  180  range-rate  measurements  over  a  20-day  tracking  period  for 
a  spacecraft  on  a  nominal  halo- type  path  near  L^.  The  standard 
deviations  of  the  three  position  states  are  reduced  to  one  measure  of 
position  uncertainty  by  taking  the  square  root  of  the  sum  of 
the  three  state  (x,  y,  z  )  variances.  The  velocity  uncertainty  is 
derived  in  a  similar  way. 

Table  2.  Comparison  of  Batch  Least  Squares  and  Kalman  Filters 


Batch  Least 

Kalman 

Percent  Deviation 

Squares (B) 

(K) 

(K-B)/K 

Range  Uncertainty  (km) 

3.96 

4.14 

4.35% 

Velocity  Uncertainty  (m/sec) 

.00220 

. 00224 

1 . 82%  1 

Clearly,  on  the  basis  of  this  limited  comparison,  both  methods  provide 
approximately  equal  estimates  of  state  uncertainty  for  equivalent  input 
errors  and  tracking  schedules. 

It  should  be  noted  that  the  results  presented  in  Table  2  are  not 
due  to  a  single  covariance  analysis  simulation.  In  fact,  all  the 
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results  presented  for  error  analysis  in  this  work  are  averages  of 
several  simulations.  For  example,  a  covariance  analysis  could  be 
completed  for  the  first  20  days  along  the  nominal  path  with  t^=  0  days; 
an.  identical  20-day  simulation  could  also  be  completed  along  the 
nominal  path  with  t^  =  10  days.  The  covariance  analysis  results  from 
these  two  runs  may  be  different  due  to  the  change  in  the  portion  of  the 
nominal  path  along  which  the  spacecraft  is  tracked,  and  average  error 
values  based  on  these  two  simulations  may  then  be  used.  The  variances 
(along  the  diagonals  of  the  covariance  matrices)  from  the  two 
simulations  for  each  of  the  states  are  averaged,  and  this  average 
variance  for  each  of  the  states  can  then  be  considered  as  the  output  cl' 
the  orbit  determination  error  analysis. 

Because  the  covariance  analysis  results  do  vary  depending  on  the 
spacecraft’s  position  in  the  orbit  (and  the  geometry  of  the  tracking 
solution),  the  covariance  analysis  results  presented  as  part  of  this 
effort  are  actually  computed  from  the  averages  of  30  simulation  runs. 
Each  run  generally  uses  a  20-day  tracking  period  unless  otherwise 
specified,  but  the  initial  epoch  for  each  simulation  is  adjusted 
(stepped  further  along  the  nominal  path).  The  adjustments  in  the 
initial  epoch  are  computed  so  that,  after  30  simulation  runs  are 
completed,  the  initial  epoch  used  on  the  last  run  has  advanced  through 
approximately  180  days  (about  one  "revolution")  of  the  orbit.  The 
positive  square  roots  (the  standard  deviations)  of  these  "average" 
variances  for  the  residual  states  can  then  be  used  to  represent  the 
state  error  levels  that  would  eventually  be  used  in  a  station-keeping 
simulation. 

An  interesting  area  for  future  research  could  be  the  coupling  of 
the  orbit  determination  error  analysis  and  station-keeping  as  one 
concurrent  simulation;  that  is,  the  error  analysis  could  produce  a 
level  of  residual  state  uncertainty  for  the  first,  say,  20  days  of 
tracking,  and  the  exact  (not  average)  uncertainty  could  then  be  input 
in  the  Monte  Carlo  simulation  of  the  station-keeping  algorithm.  The 
residual  state  uncertainty  levels  derived  from  the  next  20-day  tracking 
period  could  then,  in  turn,  be  input  at  the  40-day  time  point  as  the 
Monte  Carlo  simulation  continued.  In  general,  orbit  determination 
error  analysis  and  station-keeping  simulations  are  completed 
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separately;  this  coupling  of  the  simulations  could  produce  more 
accurate  results  in  future  research. 

The  next  step  in  this  investigation  is  the  addition  of  parameter 
errors  in  the  error  analysis.  With  the  addition  of  parameter 
uncertainty,  the  comparison  of  orbit  determination  error  analysis 
results  for  both  halo-type  and  Lissajous  orbits  can  be  completed  using 
a  larger  error  spectrum.  When  some  parameter  errors  are  assumed  to 
have  a  nonzero  level  of  uncertainty,  the  resulting  state  uncertainty 
magnitudes  for  identical  tracking  schemes  may  increase  compared  to  the 
results  of  the  previous  two  methods  (the  Kalman  and  batch  weighted 
least  squares  filter  covariance  analyses);  consider  covariance  analysis 
is  a  natural  method  to  use  when  parameter  uncertainty  is  added  to  this 
error  analysis  investigation. 


D.  Consider  Covariance  Analysis  Results 

In  the  consider  covariance  analysis  studies  completed  for  this 
work,  the  initial  purpose  was  to  attempt  to  approximately  match  the 
state  uncertainty  levels  derived  in  other  works.  The  agreement  was 
predicted  to  be  approximate  because  orbital  shapes  and  sizes  and 
dynamic  models  are  not  exactly  matched.  Also,  tracking  schedules  may 
not  be  accurately  described  for  every  study.  For  the  research  results 
compared  in  Table  3,  the  tracking  schedules  and  the  tracking  and 
parameter  uncertainty  levels  were  matched  as  closely  as  possible.  The 
position  uncertainty  is  computed  by  taking  the  square  root  of  the  sum 
of  the  first  three  elements  along  the  diagonal  of  the  aposteriori 
covariance  matrix  In  other  words,  the  total  position  uncertainty 
is  computed  by  summing  the  variances  of  the  x,  y,  and  z  states  and  then 
taking  the  square  root  of  this  sum.  The  total  velocity  uncertainty  is 
computed  in  a  similar  way. 

The  results  in  Table  3  are  arranged  in  three  groups,  separated  by 
dotted  lines.  Within  each  group,  similar  tracking  schedules  and  error 
levels  are  used  when  known;  the  covariances  computed  in  this  effort  for 
presentation  in  Table  3  are  average  values  from  30  simulations.  (The 
input  error  levels  for  several  error  analysis  studies  were  listed  in 
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Table  1.)  Eor  instance,  the  first  group  includes  works  by 

Mistretta^  and  Heuberger*  that  pertain  to  ISEE-3  studies.  Tht' 
study  by  Mistretta  is  an  error  analysis  investigation  that  simulates 
one  range  measurement  per  site  per  day  for  14  days  using  three  tracking 
sites.  It  also  assumes  one  range-rate  measurement  per  10  minutes  for 
the  2  hours  that  each  of  the  three  sites  can  adequately  view  the 
spacecraft.  An  identical  tracking  schedule  is  used  in  the  consider 
covariance  analysis  used  in  this  work  to  attempt  to  match  Mistretta’ s 
findings.  The  study  by  Heuberger  is  an  analysis  of  station-keeping 
costs  for  ISEE-3,  and  it  uses  the  error  levels  represented  in  Table  3. 

In  the  second  and  third  groupings  for  the  comparisons  of  Table  3, 
an  attempt  is  also  made  to  match  the  tracking  schedules  and  input  error 

143] 

levels  for  which  the  studies  were  conducted.  In  the  study  by  Joyce 
for  the  ISEE-3  mission  analysis,  a  tracking  schedule  of  three  range  and 
range-rate  measurements  per  site  per  day  for  each  of  three  tracking 
sites  was  simulated  for  21  days  of  tracking.  In  the  study  by 

[45] 

Rodriquez-Canabal  ,  the  tracking  schedule  includes  one  range 
measurement  per  site  per  day  and  8  range-rate  measurements  per  site  per 
day  using  two  tracking  sites  for  20  days  of  tracking. 

The  error  levels  listed  in  Table  3  show  results  that  are 
considered,  when  allowing  for  possible  force  model  and  nominal  path 
differences,  to  be  in  general  agreement  within  groupings;  thus,  with  a 
significant  level  of  confidence,  further  simulations  are  conducted  to 
quantify  the  state  uncertainty  levels  to  be  used  in  the  station-keeping 
simulations  of  chapter  three.  (An  interesting  focus  of  future  research 
might  be  to  evaluate  the  differences  in  station-keeping  costs,  if  any, 
as  a  result  of  using  the  various  orbit  determination  error  levels 
listed  in  Table  3.  This  follow-on  research  could  then  also  use  a 
statistical  "design  of  experiments"  approach  to  quantify  both  the 
sensitivity  of  the  residual  state  error  levels  calculated  in  the  error 
analysis  and,  in  turn,  of  the  station-keeping  costs  to  varying  error 
analysis  input  uncertainty  levels. 
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Table  3.  Comparison  of  Consider  Covariance  Error  Analysis  Results. 


Study 

Total  Position 
Uncertainty (km) 

Total  Velocity 
Uncertainty ( cm/sec ) 

Mistretta.*^^^ 

19.81 

1.24 

This  Work 

31.25 

3.43 

Heuberger 

42.53 

3.67 

Joyce 

s  10.0 

i  2.00 

This  Work 

14.98 

1.69 

Rodriquez- 

Canabal'*^* 

5.83 

0.44 

This  Work 

15.04 

1.69 

E.  Computing  Error  Levels  for  Station-Keeping  Simulations 

Three  methods  (the  Kalman  filter,  batch  least  squares,  and 
consider  covariance  analysis)  have  now  been  introduced  in  this  work  and 
can  be  used  to  complete  an  orbit  determination  error  analysis.  Of 
course,  other  methods  of  error  analysis  are  also  available,  but  any  of 
these  three  previously-introduced  methods  should  prove  to  be  useful. 
The  choice  between  them  is  generally  dependent  on  the  types  of  errors 
that  need  to  be  included.  (It  is  also  partly  arbitrary.)  For 
instance,  the  Kalman  filter  and  batch  least  squares  generate  basically 
equivalent  covariance  analysis  results,  given  identical  inputs  such  as 
tracking  and  initial  position  error  levels.  The  Kalman  filter  has  the 
additional  benefit  of  computing  current  covariance  matrices  and 
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position  estimates  after  every  moasui  emont  update;  batoli  least  squau'. 
depending  on  the  specific  formulation  of  the  aJgoritlim,  can 

(42) 

unfortunately  exhibit  numerical  instability  and  does  not  produce 
state  estimates  after  every  measurement  update.  Alternatively, 
consider  covariance  analysis  is  a  weighted  batch  least  squares  approach 
that  also  incorporates  parameter  uncertainty;  its  computation  is  more 
complicated,  and  the  contribution  of  parameter  uncertainty  tends  to 
accumulate  with  each  measurement  update.  (See  equation  (2-87).)  In 

[43] 

one  study  that  discussed  error  analysis  for  the  ISEE-3  mission, 
parameter  uncertainty  (specifically  solar  reflectivity  uncertainty) 
became  the  dominant  error  source  after  just  2  days  of  tracking  and,  in 
fact,  caused  the  state  uncertainty  to  increase  as  tracking  data 
continued  to  arrive.  Similar  behavior  was  noted  in  this  work  when 
solar  reflectivity  uncertainty  was  used  as  a  parameter  error  source  in 
the  consider  covariance  analysis. 

Consider  covariance  analysis  does  allow  the  inclusion  of  a  wide 
spectrum  of  uncertain  parameters  in  the  batch  weighted  least  squares 
covariance  formulation.  When  this  parameter  uncertainty  is  included  in 
the  covariance  analysis,  the  computed  residual  state  error  levels 
increase  compared  with  the  levels  that  result  if,  for  instance,  the 
Kalman  filter  ^derived  earlier)  would  be  used  in  the  error  analysis. 
However,  the  residual  state  errors  resulting  from  any  error  analysis 
study  can  be  used  in  station-keeping  simulations.  The  tracking  updates 
that  are  incorporated  in  the  Monte  Carlo  simulation  of  the 
station-keeping  algorithm  can,  in  turn,  include  the  residual  state 
uncertainty  as  one  error  source. 

The  si7.e  of  the  residual  state  errors  should  have  some  effect  on 
the  station-keeping  cost  for  a  spacecraft  in  a  libration  point  orbil. 
Several  error  sources  (such  as  planetary  mass,  solar  reflectivity,  and 
tracking  station  position  uncertainties)  have  been  shown  in  other 
studies  to  be  significant  contributors  to  residual  state  uncertainty, 
and  they  should,  therefore,  be  incorporated  in  the  evaluation  of 
station-keeping  costs.  There  are,  however,  several  ways  to  incorporate 
these  various  errors;  consider  covariance  analysis  is  not  the  only 
available  approach.  The  parameter  errors  can  be  modeled  as  random 
errors  in  a  Monte  Carlo  simulation  of  the  station-keeping  algorithm. 
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This  would  allow  the  orbit  determination  error  analysis  simulations  to 
be  completed  for  this  effort  using  only  initial  position  and  tracking 
errors  in  the  batch  least  squares-  or  Kalman  filter  covariance  analysis. 
The  residual  state  error  levels  that  are  computed  in  the  covariance 
analysis  would  undoubtably  then  be  smaller  as  the  spectrum  of  parameter 
errors  to  be  considered  is  reduced.  The  method  that  is  eventualy  used 
for  inclusion  of  parameter  errors  somewhere  in  the  station-keeping 
problem  is,  in  general,  dependent  on  the  nature  of  the  problem  and  the 
desires  of  the  investigator;  however,  the  method  used  to  model 
parameter  errors  (as  part  of  the  orbit  determination  error  analysis 
versus  their  incorporation  as  part  of  the  Monte  Carlo  simulations  of 
the  station-keeping  algorithm)  may  have  some  consistent  impact  on 
station-keeping  costs  and  could,  therefore,  be  a  fruitful  area  for 
future  research. 

For  the  current  effort,  consider  covariance  analysis  is  selected 
as  the  method  to  be  used  to  determine  the  state  uncertainty  levels  that 
will,  in  future  research,  be  used  in  station-keeping  studies.  The 
input  error  levels  (except  solar  reflectivity  uncertainty),  listed  in 
Table  1  for  this  work,  are  used'  here.  Solar  reflectivity  uncertainty 
has  been  used  in  many  of  the  previous  error  analysis  simulations 
presented  in  this  effort;  however,  several  studies^^^’^^’^°^  pertaining 
to  future  missions  involving  libration  point  orbits  do  not  use  solar 
reflectivity  uncertainty  as  an  input  to  an  orbit  determination  error 
analysis.  In  fact,  these  studies  use  the  solar  reflectivity 
uncertainty  as  a  separate  random  error  source  in  the  Monte  Carlo 
simulations  of  the  station-keeping  algorithm.  (The  station-keeping 
simulations  would  naturally  follow  completion  of  the  orbit 
determination  error  analysis  studies. ) 

In  order  to  more  closely  match  these  recent  studies,  solar 
reflectivity  uncertainty  in  this  work  is  also  input  as  an  error  source 
in  only  the  station-keeping  simulations;  however,  there  is  also  another 
reason  that  solar  reflectivity  uncertainty  is  preferred  here  as  an 
error  source  in  the  station-keeping  simulations  rather  than  in  the 
orbit  determination  error  analysis.  The  added  parameter  uncertainty, 
especially  solar  reflectivity,  causes  the  residual  state  error  levels 
to  vary  somewhat  depending  on  the  spacecraft’s  position  in  the  orbit. 


76 


The  30  covariance  analysis  trials,  completed  along  the  orbit,  can  still 
be  used  to  compute  average  residual  state  errors;  however,  the  average 
values  will  probably  not  be  as  "close"  to  the  true  residual  state  error 
values  for  a  given  20-day  segment  of  the  orbit  if  these  values  vary 
widely  along  the  orbit. 

Table  4  depicts  the  range  (lowest,  highest)  of  values  for  the 
residual  state  standard  deviations  resulting  from  a  series  of  consider 
covariance  analyses.  In  order  to  construct  this  table,  a  specific 
value  (10%,  5%,  2.5%,  or  0%)  was  selected  for  the  solar  reflectivity 
uncertainty,  and  a  series  of  30  consider  covariance  analyses  using 
20-day  tracking  periods  were  completed.  The  other  parameter 
uncertainty  levels  and  the  tracking  schedules  are  consistent  throughout 
the  simulations  used  to  construct  the  data  for  Table  4.  The  first 
covariance  analysis  is  completed  at  the  "beginning"  of  the  nominal 
halo-type  orbit  depicted  in  Figure  1-5.  (Note  that  Figure  1-5  includes 
depictions  using  only  the  three  position  states,  and  the 
"direction-of -orbit"  arrows  in  this  figure  are  positioned  at 
approximately  the  orbit  start  point.  The  state  vector  denoting  the 
"beginning"  of  the  orbit  includes  the  three  position  states  and  the 
three  velocity  states.)  Each  subsequent  covariance  analysis  is 
initiated  at  a  position  that  is  approximately  6  days  further  along 
track.  The  halo-type  nominal  orbit  is  not  perfectly  periodic,  but  it 
does  approximately  complete  a  full  revolution  in  about  180  days. 
Therefore,  this  scheme  of  error  analysis  simulations  will  cover 
approximately  one  revolution  of  the  nominal  halo-type  path. 
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Table  4.  Error  Analysis  Results  for  a  Halo-Type  Orbit:  Range  in  the 
Standard  Devia^on  for  Various  Levels  of  Solar  Reflectivity 
Uncertainty. 


Coordinate 

Minimum,. 

O'  =  10’/. 
k 

Maximum  One  Standard  Deviation  Levels 

O'  =  5%  O'  =  2.5%  ir;  =  0% 

k  k  k 

X 

(km) 

5.24, 

27.69 

4.20, 

17.96 

3.32, 

9.16 

1..38, 

1.98 

y 

(km) 

3.42, 

50.94 

2.39, 

39.43 

2.21, 

17.52 

1.68, 

3.13 

z 

(km) 

1.17, 

26.87 

1.17, 

19.03 

1.17, 

5.69 

1.16, 

5.18 

X 

(cm/sec) 

.51, 

2.55 

.51, 

1.74 

.32, 

.81 

.10, 

.36 

« 

y 

(cm/sec) 

.14, 

4.44 

.12, 

3.59 

.09, 

.98 

.08, 

.22 

• 

z 

(cm/sec) 

.12, 

1.60 

.12, 

1.25 

.11, 

.60 

.11, 

.28 

The  contributions  of  the  remaining  error  sources,  other  than  solar 
reflectivity  uncertainty,  were  not  as  noticeably  dependent  on  orbital 
position.  Consequently,  solar  reflectivity  uncertainty  is  used  as  a 
separate  error  source  in  the  station-keeping  simulations  discussed  in 
chapter  three;  it  is  not  used  to  compute  residual  state  uncertainty 
levels  that  are  the  outputs  of  an  error  analysis.  The  error  analysis, 
completed  for  later  use  in  the  station-keeping  simulations,  assumes  a 
20-day  tracking  arc  using  3  passes  per  day  from  3  separate  tracking 

[43] 

sites.  These  assumptions  closely  match  those  of  Joyce.  The  input 

error  levels  used  here  were  accurately  summarized  in  Table  1.  Mass 

parameter  uncertainty  will  also  match  the  levels  found  in  Joyce.  The 

station  position  and  tracking  uncertainty  levels  will  match  those  used 

[45] 

in  Rodriques-Canabal .  Solar  reflectivity  uncertainty  can  be  used 
in  the  station-keeping  simulations  as  a  random  variable.  The  results 
of  the  error  analysis  simulations  are  summarized  in  Table  5.  These 
error  levels  can  be  used  in  the  future  station-keeping  simulations. 
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Table  5.  Error  Levels  Produced  from  Error  Analysis  Studies. 


One  Standard 

Deviation  Levels 

Coordinate 

Halo-Type  Orbit 

Lissajous  Orbit 

X  (km) 

1.46 

1.25 

y  (km) 

2.64 

3.35 

z  (km) 

4.81 

3.19 

« 

X  (mm/sec) 

1.40 

1.25 

y  (mm/sec) 

1.85 

1.41 

2  (mm/sec) 

2.49 

2.51 

Certainly,  It  is  of  great  interest  to  compare  these  error  levels 
with  the  results  of  other  investigations  involving  spacecraft  in  halo 
(or  at.  least  halo-type)  orbits  near  the  interior  Sun-Earth  libration 
point.  Table  6  lists  the  results  of  four  such  error  analysis  studies 
that  do  not  include  solar  reflectivity  as  an  error  source. 


Table  6,  Comparison  of  Error  Analysis  Results  from  Several  Sources. 


One  Standard 

Deviation  Error  Levels 

Coordinate 

r  45  ] 

Rodriquez-Canabal 

Simo'*’^ 

c..  -ISO] 

Simo 

This  Work 

X  (km) 

2.7 

1.5 

1.73 

1.46 

y  (km) 

3.9 

2.5 

2.24 

2.64 

z  (km) 

3.4 

15.0 

5.48 

4.81 

X  (mm/sec) 

2.4 

1.0 

1.41 

1.40 

y  (mm/sec) 

3.5 

1.0 

1.41 

1.85 

z  (mm/sec) 

1.3 

3.0 

2.45 

2.49 
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The  error  analysis  results,  computed  in  this  work  for  use  in  the 
station-keeping  simulations,  agree  most  closely  with  the  levels 
computed  in  Simo  .  The  agreement  can  be  clearly  seen  in  Table  6. 
(It  should  be  noted  that  there  are,  of  course,  small  differences  in  the 
nominal  paths  and  force  models  used  for  the  works  listed  in  Table  6. ) 

When  all  the  various  input  error  sources  are  combined  over  the  20 
days  of  tracking,  the  error  analysis  predicts  the  levels  of  state 
uncertainties  listed  in  Table  5.  Technically,  these  error  magnitudes 
do  appear  to  differ  for  halo  and  Lissajous  orbits;  however,  two 
different  Lissajous  nominal  orbits  may  show  the  same  variations. 
Station-keeping  costs,  determined  through  simulations  using  the  above 
error  levels,  may  or  may  not  differ.  The  testing  of  the  significance 
of  these  differences  in  the  error  levels  will  be  undertaken  in  future 
station-keeping  simulations;  however,  these  Monte  Carlo  simulations 
will  require  assumptions  about  both  the  type  of  probability 
distribution  (Gaussian  or  uniform,  for  instance)  and  the  mean  value  of 
the  errors.  Therefore,  one  last  area  of  inquiry  in  this  error  analysis 
study  should  be  concerned  with  the  shape  (distribution)  and  position 
(means)  of  the  probability  distributions  for  the  random  errors  listed 
in  Table  5. 


F.  Probability  Distribution  of  the  Resulting  Errors 

Most,  if  not  all,  station-keeping  studies  use  random  errors  as  if 
they  were  Gaussian  with  a  zero  mean.  That  is,  the  Monte  Carlo 
simulations  associated  with  the  control  algorithms  use  standard 
deviation  levels  such  as  those  in  Table  5  in  a  random  number  generator 
for  which  the  probability  distributions  and  the  mean  values  of  the 
errors  must  be  assumed.  Generally,  the  uncertainty  is  assumed  to  be 
"white  noise" — errors  that  are  zero  mean  and  that  tend  to  closely 
follow  the  normal  distribution.  A  covariance  analysis  is  generally  a 
linear  algorithm  and,  in  a  linear  analysis,  sums  of  Gaussian  random 
variables  are  clearly  also  Gaussian.  In  this  case,  a  linear  error 
analysis  that  then  assumes  Gaussian,  zero  mean  input  errors  will,  in 
turn,  compute  residual  state  errors  that  are  also  zero  mean  and 
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Gaussian.  However,  it  is  not  so  clear  what  distribution  the  residual 
state  errors  might  follow  if  the  random  errors  were  alternatively 
propagated  using  the  nonlinear  equations  of  motion,  a  method  sometimes 
used  in  this  work. 

For  this  portion  of  the  research,  the  input  error  levels  listed 
for  this  work  in  Table  1  are  assumed,  and  the  nonlinear  equations  of 
motion  are  integrated  forward  for  20-day  arcs.  Range  and  range-rate 
tracking,  station  position,  mass  parameter,  and  solar  radiation 
pressure  random  errors  were  simulated  using  zero-mean, 
normally-distributed  variables.  Numerical  integration  routines  and 
random  number  generators  available  in  the  software  package 
386-Matlab  were  used  in  this  effort.  Three  hundred  random  trials 
were  conducted.  ilistograms  of  the  resulting  errors  in  position  and 
velocity  relative  to  the  nominal  orbit  for  each  of  the  six  states  are 
displayed  in  Figures  3-5  through  3-10.  The  histograms  use  ten  classes 
(of  equal  class  width)  with  midpoints  plotted  along  the  horizontal 
axis. 

The  vertical  axis  represents  the  frequency  of  responses  counted  in 
the  classes.  These  six  depictions  show  distributions  of  errors  that 
appear  to  be  normally  distributed  with  a  mean  of  approximately  zero. 
(The  appearance  of  a  histogram  associated  with  an  exactly  normally 
distributed  random  variable  will  be  discussed  at  the  end  of  this 
chapter. )  Statistical  tests  can  be  conducted  to  verify  these 
assumptions  of  the  distribution  type  and  location  of  the  mean  value. 

Initially,  it  must  be  determined  if  the  distributions  are,  indeed, 
Gaussian.  Thus,  the  hypothesis  that  the  distributions  are  normal 
(assumed  true)  versus  the  hypothesis  that  they  are  not  is  to  be  tested. 

(51  j 

It  must  be  proven  to  a  confidence  level  of,  say,  95%.  Here  it  is 

2 

appropriate  to  use  chi-squared  (;\;  )  goodness  of  fit  tests  to  tost  the 
hypotheses.  The  goodness  of  fit  test  is  conducted  separately  for  each 
of  the  six  coordinates.  It  compares  the  observed  frequencies  (f^)  in 
the  classes  used  to  construct  the  histograms  shown  in  Figures  3-5 
through  3-10  to  the  expected  frequencies  (F^)  under  the  hypothesis  that 
the  distribution  is  normal. 
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Figure  3-5.  Histogram  of  the  x  Excursion  From  the  Integrated  Path. 
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y  excursions  in  kilometers  from  the  integrated  path 


Figure  3-6.  Histogram  of  the  y  Excursion  From  the  Integrated  Path. 
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Figure  3-7.  Histogram  of  the  z  Excursions  From  the  Integrated  Path 
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Figure  3-9.  Histogram  of  the  y  Velocity  Excursion  From  the  Integrated  Path. 
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Figure  3-10.  Histogram  of  the  z  Velocity  Excursion  From  the  Integrated  Path. 
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2 

For  these  tests,  the  test  statistic  X  is  defined 


1=1 


where  k  denotes  the  number  of  classes  in  the  histogram  for  each  of  the 

six  variables.  Here,  k  is  ten  (or  nine  If  classes  must  be  pooled),  and 

the  degrees  of  freedom  for  the  goodness  of  fit  test  are  derived  from 

the  number  of  classes.  In  this  investigation,  the  degrees  of  freedom 

for  the  chi-squared  test  will  be  reduced  by  one  because  the  total  of 

the  expected  frequencies  (f^)  must  be  equal  to  the  sample  size.  It  is 

further  reduced  by  two  because  both  the  mean  and  variance  used  to 

compute  the  expected  frequencies  are  computed  from  the  sample  data. 

The  X  random  variable  for  these  tests  then  has  k-2-l=7  (or  6  if 

pooled)  degrees  of  freedom,  a  95%  level  of  confidence,  and  is  equal  to 

14.07  (or  12.59  if  pooling  is  necessary).  This  value  for  x  is  found 

[511 

in  a  table  of  values  for  the  distribution. 

The  hypotheses  and  the  decision  rule  used  for  each  of  the  six 
tests  are: 


Hypotheses: 

H^:  The  Probability  Distribution  for  (x,y,z,x,y,  or  z)  is  normal. 
H^:  The  Probability  Distribution  is  not  normal. 

Decision  Rule: 

If  s  x^>  conclude  H^. 

If  i  x^>  conclude 

The  results  of  these  tests  are  summarized  in  Table  7  for  the  Monte 
Carlo  simulation  of  300  runs. 
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Table  7.  Results  of  the  Error  Analysis  Goodness  of  Fit  Tests. 


Error 

2 

X 

Conclusion 

X 

1.712 

12.59 

H  : Gaussian 
0 

y 

12.486 

14.07 

H  ; Gaussian 
0 

z 

8.524 

14.07 

H  : Gaussian 
0 

X 

5.430 

12.59 

H  : Gaussian 
0 

• 

y 

12.338 

12.59 

H  : Gaussian 
0 

• 

z 

9.725 

12.59 

i)  ;  Gaussian 
0 

V 


r 


The  general  conclusion  of  this  hypothesis  test  is  that  the 
residual  state  errors  that  result  from  the  error  analysis  can  be 
appropriately  modeled  as  Gaussian.  Hence,  the  state  errors  used  in  the 
Monte  Carlo  simulations  of  the  station-keeping  algorithms  can  be 
modeled  as  normally  distributed  random  variables.  (It  should  be  noted 
that  the  results  of  the  chi-squared  goodness  of  fit  test  can  be 
presented  in  different  ways,  and,  in  fact,  other  types  of  goodness  of 
fit  tests  are  also  available.  For  instance,  the  well-known  Statistical 
Analysis  Software  from  the  SAS  Institute  uses  a  Shapiro/Wilk  test  for 
normality  and,  when  the  same  residual  state  error  data  is  used  as 
inputs  to  this  software,  the  program  indicates  that  each  of  the  errors 
is  strongly  Gaussian  with  a  mean  of  zero.) 

It  must  also  be  determined  if  these  normally  distributed  state 

errors  are  also  zero-mean.  Now  the  hypothesis  to  be  tested  is  that  the 

population  means  for  the  random  errors  are  zero.  This  hypothesis  is 

tested  six  independent  times — once  for  each  state  variable.  The 

population  means  are  denoted  as  u  ,  u  ,  u  ,  u*  ,u-  ,ii*,  and  the  conduct 

X  y  z  X  y  z 

of  the  hypothesis  tests  assumes  that  these  population  means  are  zero 

unless  proven  otherwise.  For  this  statistical  test,  the  standard 

normal  distribution  (Z)  is  the  assumed  working  distribution.  The 

statistical  test  statistic  is  denoted  by  Z  and  is  calculated  by  (using 
the  residual  state  x  as  a  specific  example): 
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sample  mean 


X 


X 


Z 


standard  error  of  the  mean 


s{x} 


s{x}/(n'^^) 


where  s{x}  is  the  standard  deviation  of  the  fandom  variable  in  question 
(which  is  X  here)  and  n  is  the  random  sample  size.  Here,  n  is  300 
random  trials.  The  level  of  confidence  used  for  this  test  is  again 
9554,  but  here  the  risk  (5%)  is  divided  into  both  tails  of  the 
distribution.  Therefore,  Z  =  1,960,  which  can  be  found  in  a  table  of 
values  for  the  standard  normal  distribution. 

The  alternatives  ana  decision  rule  are  then: 


Hypotheses: 

=  0. 

0  X 

H  :p  ^0. 

l  x 

Decision  Rule: 

» 

If  |Z  I  a  Z  ,  conclude  H^. 

If  |Z*|  a  Z  ,  conclude  H^. 

The  results  are  summarized  below  in  Table  8. 


Table  8.  Hypothesis  Test  for  Zero  Error  Means. 


» 

Error 

-Z 

Z 

Z 

Conclusion 

X 

-1.960 

-1.051 

1.960 

H  :m  =  0. 

0  X 

y 

-1.960 

1.193 

1.960 

H  :m  =  0. 

0 

z 

-1.960 

-.403 

1.960 

H  :p  =  0. 

0  2 

X 

-1.960 

-.679 

1.960 

H„:p-=  0. 

0  X 

y 

-1.960 

1.084 

1.960 

H  :p-=  0. 

0 

z 

-1.960 

-.518 

1.960 

H  :p-=  0. 

0 
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It  may  be  somewhat  Interesting  to  discuss  the  differences  between 
a  distribution  that  appears  to  be  Gaussian  and  one  that  is,  in  fact, 
Gaussian.  The  most  appropriate  way  to  present  this  comparison  is 
through  two  histograms,  one  for  a  variable  tested  in  the  discussion 
above  and  one  for  a  Gaussian  distribution  with  an  identical  mean  and 
standard  deviation.  (The  shape  of  a  Gaussian  distribution  is 
determined  by  its  standard  deviation;  the  location  of  the  distribution 
on  the  real  number  line  is  determined  by  the  mean. )  The  residual  y 
state  variable  in  both  the  goodness  of  fit  test  and  the  hypothesis  test 
for  a  mean  of  zero  provided  the  results  that  were  "closest  to  causing 
rejection"  and  could  thus  be  considered  to  depart  most  from  the 
Gaussian  distribution.  The  mean  of  y  is  .  156635  and  the  standard 
deviation  of  y  is  3.34769  in  both  histograms  of  Figures  3-11  and  3-12. 
The  histogram  in  Figure  3-11  is  for  the  actual  residual  y  state,  and 
the  observed  frequencies  f^  are  included  in  each  class;  the  histogram 
in  Figure  3-12  is  of  a  truely  Gaussian  distribution,  and  the  expected 
frequencies  F^  are  included  in  each  class. 

With  the  completion  of  both  types  of  statistical  tests  for  the 
residual  state  errors,  the  state  errors  can  be  treated  as  zero-mean 
Gaussian  random  variables.  Follow-on  research  incorporating 
station-keeping  simulations  may  use  these  characteristics  to  model  the 
random  errors,  and  the  magnitudes  of  the  residual  state  errors  used  for 
the  Lissajous  and  halo-type  orbit  simulations  may  be  as  depicted  in 
Table  5.  The  solar  reflectivity  uncertainty  can  be  input  as  a  separate 
random  variable  with  a  standard  deviation  of  13%,  as  depicted  in  Table 
1.  In  station-keeping  simulations,  the  propellant  used  for  a  period  of 
station-keeping  can  be  treated  as  a  random  variable  and  subjected  to 
statistical  tests  similar  to  those  used  in  this  section. 
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y  excursions  in  kilometers  from  the  integrated  path 


mean  =  .1566 

standard  deviation  =  3.3477 


Figure  3-1 1 .  Histogram  of  y  Excursion  From  Integrated  Path  Including  Observed  Frequencies. 
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mean  =  .1566 

standard  deviation  =  3.3477 


Figure  3-12.  Histogram  of  y  Excursion  From  the  Integrated  Path  Using  Gaussian  Expectations. 
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